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Abstract 


Gouging  is  a  type  of  structural  failure  that  becomes  important  when  two  metals 
slide  against  each  other  at  velocities  in  the  range  of  1.5  kilometers  per  second.  A 
computer  model  has  previously  been  used  to  model  the  development  of  gouging  at  the 
Holloman  High  Speed  Test  Track.  This  model  has  not  been  experimentally  verified  to  be 
correct,  due  to  the  complexity  of  the  model.  This  research  develops  a  simplified  model 
that  can  be  experimentally  verified. 

The  computer  program  utilized  in  this  research  was  studied  to  determine  the  most 
appropriate  options  to  use  in  simulations.  This  was  accomplished  by  modeling  a  Taylor 
impact  test  and  comparing  to  published  experimental  results. 

The  cylindrical  impact  specimen  utilized  in  the  simplified  model  was  developed 
through  use  of  the  Buckingham-Pi  theorem,  and  can  be  fired  from  most  standard 
compressed  air  guns.  Simulations  using  the  simplified  model  showed  excellent 
agreement  with  simulations  using  the  physical  sled  properties.  Plasticity  observed  in 
both  the  rod  and  target  was  very  similar  to  that  seen  in  the  physical  sled  simulations.  The 
high-  pressure  core,  which  initiates  gouging  in  the  physical  sled  simulation,  was  found  to 
exist  in  the  simplified  model  as  well. 
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AN  INVESTIGATION  OF  A  SIMPLIFIED  GOUGING  MODEL 


Chapter  1  -  Introduction 

In  recent  years,  the  United  States  Air  Force  (USAF)  has  spent  considerable 
money  and  effort  investigating  methods  to  improve  the  understanding  of  very  high 
velocity  impact  testing.  The  Air  Force  utilizes  a  sophisticated  test  facility  located  at 
Holloman  Air  Force  Base  (AFB)  in  New  Mexico  to  perform  the  majority  of  the  advanced 
tests  that  they  perform.  The  main  system  consists  of  a  narrow  gauge  rail  system 
approximately  ten  miles  long,  which  is  used  to  guide  a  rocket  sled,  see  Figure  1,  at 
extremely  high  velocities.  In  April  2003,  the  Holloman  High  Speed  Test  Track  (HHSTT) 
achieved  a  world  record  velocity  of  2884.9  m/s.  The  846th  Test  Squadron,  which 
operates  the  HHSTT,  is  working  to  increase  the  maximum  velocity  to  10,000  ft/s  or 
approximately  3  km/s. 


Figure  1-  Rocket  Sled  at  HHSTT 
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In  order  to  increase  the  maximum  velocity  of  the  test  track,  engineers  at  the  846th 
Test  Squadron  and  the  Air  Force  Institute  of  Technology  (AFIT)  are  investigating 
methods  that  will  mitigate  the  phenomenon  known  as  gouging.  The  rocket  sled  seen  in 
Figure  1  is  attached  to  the  narrow  gauge  rail,  which  is  made  of  1080  steel,  via  four 
slippers,  which  are  made  of  high  strength  VascoMax  300  steel,  seen  schematically  in 
Figure  2.  Gouging  typically  occurs  at  velocities  greater  than  1.5  km/s,  and  can  range  in 
severity  from  a  shallow  gouge  more  reminiscent  of  damage,  to  a  total  structural  failure  of 
the  sled/rail  assembly. 


According  to  Laird  [1],  gouging  occurs  “when  inertial  forces  are  so  great  that  the 
materials  exhibit  fluid  like  behavior.  Shock  induced  pressure  creates  a  region  of 
plasticity  under  the  location  of  impact.  Tangential  motion  of  one  body  with  respect  to  the 
other  deforms  or  shears  material  at  these  points  and  results  in  deformation  of  the  parallel 
surfaces  that  impinge  on  each  other  in  a  continuous  interaction.”  Gouges  occur  due  to  the 
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fact  that  gaps  exist  between  the  slipper  and  the  rail,  which  allows  the  rocket  sled  to  move 
perpendicular  to  the  rail  direction,  as  well  as  to  yaw,  pitch,  and  roll.  Any  of  these  relative 
motions  can  create  a  situation  whereby  the  slipper  impacts  the  rail  at  an  extremely  high 
velocity,  which  in  some  cases  can  cause  gouging  to  occur.  A  typical  gouge  is  seen  in 
Figure  3.  [2] 


Direction  of 
Sled  Travel 


Figure  3  -  Gouged  Rail  Section 

A  gouge  in  the  rail  material  is  characterized  by  a  ridge  of  material  in  front  of  the 
initiation  point,  followed  by  a  steep  drop  into  the  gouge,  which  gradually  gets  shallower 
as  the  gouge  continues  to  form,  until  the  slipper  loses  contact  with  the  rail  and  the 
gouging  is  complete.  The  schematic  view  of  a  gouge  seen  in  Figure  4  shows  the  classic 
teardrop  shape  common  to  most  gouges. [3]  A  typical  gouge  is  approximately  15.0  cm 
long  by  around  .6  cm  deep. 


Direction  of 
Sled  Travel 


PLAN  VIEW 


CROSS  SECTION 


Figure  4  -  Schematic  View  of  Gouge 
3 


The  question  of  the  magnitude  of  the  vertical  velocity  with  which  the  slipper 
strikes  the  rail  was  answered  by  Hooser  [4]  via  a  simulation  of  a  10,000  ft/s  (3048  m/s) 
rocket  sled  test  using  the  Dynamic  Analysis  and  Design  System  (DADS).  Hooser  found 
that  the  vertical  velocities  tended  to  be  on  the  order  of  1-2  m/s  when  a  horizontal  velocity 
of  1.5  km/s  was  used.  This  shows  that  the  impacts  occurring  during  a  typical  HHSTT 
experiment  are  at  angles  of  approximately  0.03  degrees. 

Both  Laird  and  Szmerekovsky  have  studied  the  slipper  and  rail  interactions 
analytically  using  the  Sandia  National  Laboratory  hydrocode  CTH.  The  Szmerekovsky 
model  is  accepted  as  the  standard  simulation  of  the  gouging  phenomena,  because  he  used 
the  actual  test  sled  conditions  to  perform  the  simulation.  The  work  of  these  authors 
proved  very  useful  in  understanding  and  modeling  the  physics  of  what  was  occurring  as  a 
gouge  forms,  however  due  to  the  complexities  of  the  models  that  were  used,  they  were 
not  able  to  be  verified  experimentally.  A  model  which  captured  the  important  gouging 
phenomena,  yet  was  simple  enough  to  be  verified  experimentally  would  be  helpful,  as  it 
would  allow  simple  simulations  and  experiments  which  would  give  operators  at  HHSTT 
an  idea  of  what  conditions  would  create  gouging.  Knowing  when  gouging  would  occur 
would  allow  operators  to  avoid  certain  test  configurations  or  slipper  geometries,  and 
possibly  allow  the  design  of  a  slipper  that  would  eliminate  gouging. 
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Chapter  2  -  Background  and  Theory 


2.1  Conservation  Equations 


Virtually  all  fields  of  mechanics  and  dynamics  are  based  upon  the  same  three 
fundamental  conservation  laws:  the  conservation  of  mass,  the  conservation  of 
momentum,  and  the  conservation  of  energy,  however  they  become  extremely  important 
when  the  problem  is  dependent  upon  inertia.  These  equations  are  presented  here  for 
completeness. 

Conservation  of  mass  states  that  mass  cannot  be  created  nor  destroyed. 
Mathematically  this  is  stated 

|  pdV  =  constant  (1) 

v 


where  p  is  density  and  V  is  the  volume  of  the  material. 

The  conservation  of  momentum  can  be  stated  in  a  number  of  different  ways. 
Perhaps  the  simplest  statement  is  that  force  equals  mass  times  acceleration,  or 


77  dv 

t  =  m  — 
dt 


(2) 


where  F  is  the  force  applied,  m  is  the  mass  acted  upon,  v  is  the  velocity  of  the  mass,  and  t 
is  the  time  over  which  the  event  occurs.  Another  useful  statement  of  the  conservation  of 
momentum  is  the  impulse-momentum  relation,  which  is  obtained  by  multiplying  both 
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sides  of  the  conservation  of  momentum  equation  by  dt  and  integrating  over  a  period  of 
time,  giving 

/  =  J  Fdt  =  J  mdv  =  mv  f  -  mva  (3) 

where  /  is  the  impulse  applied  over  some  period  of  time  by  the  applied  force,  v/and  v0  are 
the  initial  and  final  velocities  of  the  mass,  and  therefore  the  right  hand  side  of  the 
equation  is  the  momentum  change  over  some  time  period.  This  represents  the 
conservation  of  linear  momentum.  A  similar  equation  for  the  conservation  of  angular 
momentum  exists  as  well,  which  leads  to  the  fact  that  the  stress  tensor  is  symmetric, 

<j  =  a  . 

v  j' 

Lastly,  the  conservation  of  energy  equation  written  for  a  discrete  set  of  j  masses  is 


+  Wr 


0->l 


(4) 


where  E  is  the  internal  energy  source,  the 


]_ 

2 


p\r  terms  are  the  kinetic  energy,  and  W0^x  is 


the  work  done  on  the  system.  The  subscript  0  represents  the  initial  state  and  the  subscript 
1  represents  the  final  state. 


2.2  Stress  Waves  in  a  Continuum 

Stress  waves  develop  in  materials  every  time  that  a  pressure  is  applied.  In 
continuum  mechanics  this  fact  is  often  ignored  and  the  loading  is  said  to  be  either  static 
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or  quasi-static.  This  practice  works  well  in  most  structural  analysis  problems  where 
pressure  or  loads  are  applied  very  slowly,  but  in  higher  velocity  impact  dynamics 
problems,  this  assumption  is  a  poor  simplification  of  what  actually  occurs. 

A  stress  wave  develops  when  a  load  is  applied  instantaneously;  the  material 
particles  immediately  surrounding  the  application  point  initially  support  all  of  the 
pressure  via  local  acceleration.  As  the  affected  particles  accelerate,  they  build  up 
compressive  stresses  in  the  particles  further  from  the  application  point.  The  material 
particles  behind  the  wave  begin  to  accelerate  as  well,  as  the  particles  that  were  initially 
struck  continue  to  move  to  support  the  pressure  applied,  compressive  stresses  develop 
behind  the  wave.  When  the  compressive  stress  is  equivalent  to  the  applied  pressure,  the 
relative  motion  can  stop,  because  the  pressure  applied  is  supported  by  the  compressive 
stress  that  has  been  developed.  This  process  continues  further  away  from  the  impact 
point,  until  the  whole  material  experiences  the  stress  wave. [5]  The  motion  of  these  so 
called  particles  is  what  is  governed  by  the  conservation  of  momentum  equation,  which 
allows  the  calculation  of  the  particle  velocities,  which  is  dependent  upon  the  impact 
velocity.  The  speed  of  propagation  of  the  stress  wave  however  will  be  shown  to  be  a 
material  property. 

A  dynamic  impact  event  can  create  stress  waves  of  widely  varying  intensities  and 
velocities  depending  upon  the  impact  velocity.  At  very  low  velocities,  an  elastic  wave 
will  form  which  defonns  the  rod  elastically.  At  higher  velocities  an  elastic  wave  will 
form,  followed  by  a  slower  moving  plastic  wave,  which  will  plastically  deform  a  portion 
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of  the  material.  At  velocities  greater  than  the  wave  propagation  velocity  of  the  material, 
generally  2  km/s,  a  shock  wave  is  formed. [5] 


2.2.1  Elastic  Stress  Waves 

Simplifying  the  general  impact  case  to  one  of  a  uniaxial  rod  impact,  it  is  possible 
to  calculate  the  magnitude  of  the  elastic  stress  wave  applied  to  the  rod.  Assume  that  an 
initially  stationary  rod  is  impacted  by  a  rigid,  semi-infinite  plate  moving  with  a  constant 
velocity  of  v0,  which  is  less  than  the  material  sound  speed  c,  from  the  left,  as  seen  in 
Figure  5.  After  the  impact,  an  elastic  stress  wave  develops  which  travels  to  the  right  at 
the  material  sound  speed.  Behind  the  stress  wave,  the  particles  have  been  accelerated  to  a 
velocity  equal  to  the  impact  velocity  vQ. 
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Figure  5  -  Rod  Impact  Experiment 

An  impulse-momentum  balance  shows  that  all  material  behind  the  stress  wave  is 
in  motion,  and  therefore  the  momentum  of  that  material  is 
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pAcAtv() 


(5) 


where  cAt  is  the  length  of  the  rod  which  is  moving  with  a  particle  velocity  of  va  and  A  is 
the  cross  sectional  area  of  the  rod.  Therefore  Ac  At  is  the  volume  of  material  moving  with 
a  particle  velocity  of  va. 

The  impulse  is  known  to  be  the  integral  of  the  force  over  a  period  of  time.  In  this 
uniaxial  case,  the  force  is  the  stress  times  the  area  over  which  it  acts,  giving  an  impulse  of 

oAdt  (6) 

where  a  is  the  compressive  stress  occurring  due  to  the  passage  of  the  stress  wave. 
Applying  the  conservation  of  momentum  to  these  two  equations,  assuming  an 
infinitesimal  time  step,  and  dividing  both  sides  by  the  area  and  the  time  step  gives  the 
magnitude  of  the  elastic  compressive  stress  wave  as 

<7  =  P°Vo  (7) 

where  c  is  the  material  sound  speed. 

The  speed  of  sound  in  a  material  is  a  direct  result  of  the  wave  equation 
encountered  in  differential  equations.  To  begin,  imagine  an  element  of  material  through 
which  a  disturbance  is  passing,  see  Figure  6.  The  edges  of  the  element  are  at  positions  of 
x  and  x+dx  as  measured  from  a  fixed  coordinate  system. 

Assuming  that  tension  is  positive  in  the  positive  x  direction,  the  stress  on  the  left 
and  right  sides  of  the  element  are 
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-  (7 


x 


and 


cr 


(8) 


respectively.  The  conservation  of  momentum  for  the  infinitesimal  element  of  area  A  is 
given  by 


A(-  <jx  )dt  +  A 


cr,.  + 


dx 


dt  =  (pA)d' 


v 


(9) 


where  the  left  side  of  the  equation  is  the  impulse  applied  due  to  the  stress  on  the  left  and 

right  sides,  and  the  right  side  of  the  equation  is  the  momentum  imparted  over  a  time  step. 

dx 

t 


1 - > 

x 

Figure  6  -  Material  Element 
Dividing  both  sides  by  Adt  gives 


dt rv  dv 

—  =P  — . 

dx  dt 


(10) 


Utilizing  the  knowledge  that  strain  and  velocity  are 


£ 


du 

dx 


and 


(11) 
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it  can  be  proven  that 


ds  _  dv 
dt  dx ' 


(12) 


If  it  is  assumed  that  stress  is  a  function  only  of  strain,  then  the  previous  equations  can  be 
combined  to  give  the  classic  one-dimensional  wave  equation, 


2  d2u  _  d2u 

a?- a? 


(13) 


where 


da/ 

c2  (s)=^3k. 
P 


(14) 


For  the  case  where  the  impact  is  elastic,  the  numerator  is  actually  the  elastic  modulus  of 
the  material,  which  gives  the  wave  velocity  in  the  case  of  an  elastic  impact  as, 


C  = 


(15) 


where  Em  is  the  elastic  modulus  of  the  material  and  p  is  the  density  of  the  material.  [5] 

It  is  also  interesting  to  understand  how  waves  propagate  through  a  medium,  and 
what  occurs  at  material  boundaries.  This  information,  along  with  the  definition  of  stress 
in  a  one-dimensional  impact,  will  be  used  later  to  verify  an  impact  model  developed. 
Assume  that  a  homogenous  metal  rod,  Figure  7,  strikes  a  semi-infinite  rigid  wall  at  a 
velocity  of  v0  perpendicular  to  the  wall.  In  this  case  we  will  assume  that  there  are  no 
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three-dimensional  effects  due  to  the  impact,  that  the  impact  is  completely  elastic,  and  that 
the  wave  velocity  within  the  rod  is  constant.  Immediately  before  impact,  the  rod  is 
moving  at  a  uniform  velocity  of  v0  and  the  internal  stress  is  uniformly  zero. 


Length  =  L 


Figure  7  -  Rod  Impact  Experiment 

At  the  moment  of  impact,  a  stress  wave  is  formed  which  travels  to  the  right  at  the 
wave  velocity  of  the  material,  here  repeated  as 


To  the  left  of  the  wave,  a  constant  particle  stress  of 

&  =  PCV0  (17) 

is  developed  throughout  the  rod.  To  maintain  continuity  at  the  boundary  of  the  impact,  it 
is  required  that  the  materials  on  either  side  of  the  boundary  must  have  the  same  velocity. 
Therefore,  the  impacted  end  of  the  rod  must  immediately  slow  down  to  a  velocity  of  zero, 
because  the  rigid  wall  has  no  velocity,  Figure  8.  While  the  left  end  of  the  rod  has  a  zero 
velocity,  the  right  end  of  the  rod  continues  to  travel  to  the  left  at  the  initial  velocity, 
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because  the  right  end  hasn’t  “felt”  the  impact  on  the  left  end,  because  the  speed  of  sound 
in  the  material  limits  how  fast  a  disturbance  can  be  “felt”.  This  fact  will  be  used  later  in 


verifying  a  CTH  model. 


c 


r 

o! 

M  II 

>  b 

O  o 

>  II 

r 

to<t<L/c 

Figure  8  -  Rod  Impact  Experiment  before  Reflection 
At  some  later  point  in  time,  the  stress  wave  has  traveled  the  full  length  of  the  rod, 

and  has  reached  the  right  end.  Since  the  rod  has  a  known  length  of  L,  it  is  also  known 

that  the  wave  will  reach  the  right  end  at  t  =  —  .  Since  the  wave  has  traveled  the  length  of 

c 

the  rod,  the  whole  rod  is  now  under  a  constant  particle  compressive  stress  as  defined 
earlier.  Additionally,  the  rod  has  momentarily  stopped  since  the  kinetic  energy  of  impact 
has  been  turned  into  internal  strain  energy,  Figure  9. 


Figure  9  -  Rod  Impact  Experiment  at  Reflection 
From  mechanics  of  materials  it  is  known  that  a  free  surface  cannot  support  an 

applied  stress,  therefore  after  the  reflection  of  the  stress  wave,  the  right  side  must  have 
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zero  stress.  This  means  that  when  a  stress  wave  reflects  from  a  free  surface,  it  reflects 
with  the  opposite  sign  of  what  it  started  out  as.  Therefore  the  initially  compressive  stress 
wave  reflects  as  a  tensile  wave,  effectively  zeroing  the  stress  behind  the  now  left  running 
wave.  Since  the  stress  wave  was  initially  moving  to  the  right,  and  there  is  no  constraint 
on  the  right  side,  the  right  hand  side  of  the  rod  continues  to  travel  to  the  right,  while  the 
tensile  stress  wave  continues  to  the  left,  Figure  10.  This  fact  will  also  be  used  to  verify  a 
CTH  model. 
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Figure  10  -  Rod  Impact  Experiment  after  Reflection 
When  the  stress  wave  returns  to  the  left  hand  end  of  the  rod,  it  is  attempting  to 

apply  a  tensile  load  to  the  rigid  wall,  which  is  impossible  since  the  materials  are  not 

physically  joined  in  any  manner.  Since  the  tensile  stress  can’t  be  supported  by  the 

interface,  the  rod  then  separates  from  the  wall  and  rebounds  away  at  the  initial  impact 

velocity,  unstressed,  Figure  11. 
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Figure  11  -  Rod  Impact  Experiment  after  Unloading 
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An  alternative  way  to  understand  the  progression  of  particle  stress  waves  is 
through  the  use  of  a  Lagrangian  diagram,  which  tracks  the  position  of  stress  waves  within 
a  rod  as  time  progresses.  The  Lagrangian  diagram  for  the  rod  impact  problem  just 
presented  is  shown  in  Figure  12. 


Figure  12  -  Lagrangian  Diagram  for  Rod  Impact  Problem 
At  the  beginning  of  the  impact  event,  the  rod  is  said  to  be  the  x-axis  from  x=0  out  to  x=L. 

In  region  1 ,  the  rod  is  unstressed,  and  traveling  with  a  constant  particle  velocity  of  v0  to 

the  left.  At  the  moment  of  impact  at  point  “a”,  a  compressive  particle  stress  wave  of 

magnitude  <J  —  pcvo  is  formed,  and  travels  to  the  right  through  region  2.  As  was 

described  earlier,  the  particle  velocity  at  the  impacted  end  must  be  equal  to  zero,  because 
the  wall  is  assumed  infinite;  therefore  in  region  2,  the  particle  velocity  is  zero.  When  the 
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particle  stress  wave  reaches  the  right  side  of  the  rod  at  point  “b”,  it  must  reflect  as  a 
tensile  wave  of  magnitude  <J  =  pcvo  ,  because  it  must  leave  the  free  edge  on  the  right 

stress  free.  The  compressive  wave  moving  to  the  right  imparts  a  particle  velocity  of  v0  to 
the  right.  When  the  tensile  wave  reaches  the  wall  at  point  “c”,  it  attempts  to  apply 
tension  to  the  wall,  but  finds  that  there  is  no  resistance  to  this  tension,  and  therefore  the 
rod  pulls  away  from  the  wall,  and  moves  off  to  the  right,  in  an  unstressed  state. 

2.2.2  Plastic  Stress  Waves 

The  stress-strain  curve  for  most  materials  is  characterized  by  a  linear  elastic 
portion  at  low  strains,  followed  by  another  region  that  may  or  may  not  be  linear  as  well, 
see  Figure  13.  The  point  where  stress-strain  behavior  stops  being  strictly  elastic  is  known 
as  the  elastic  limit,  dynamic  yield  stress,  or  Hugoniot  elastic  limit.  When  a  material  is 
acted  upon  by  a  stress  greater  than  this  limiting  value  of  stress,  the  material  plastically 
deforms.  In  an  impact  experiment,  the  initial  stress  wave  is  an  elastic  stress  wave,  which 
is  followed  by  a  plastic  stress  wave  that  initiates  when  the  particle  stress  applied  reaches 
the  dynamic  yield  stress.  This  plastic  stress  wave  is  slower  than  the  elastic  stress  wave. 
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Figure  13  -  Generic  Stress-Strain  Curve 
Two  basic  theories  have  been  developed  for  describing  the  propagation  of  a 

plastic  wave.  The  theories  differ  in  their  assumption  of  the  importance  of  the  strain  rate 

in  determining  material  response.  The  first  theory  for  describing  the  behavior  of 

materials  undergoing  plastic  deformation  is  called  the  rate-independent  theory.  As  the 

name  implies,  this  theory  makes  the  assumption  that  there  is  a  single  dynamic  stress- 

strain  curve  that  describes  material  behavior,  regardless  of  the  rate  at  which  the  strain  is 

applied.  This  theory  arose  from  the  assumption  that  a  material  had  a  bilinear  stress-strain 

curve,  as  seen  in  Figure  14.  This  theory  predicted  two  distinct  wave  fronts  would  be 

formed,  and  would  propagate  through  the  material  at  distinct  velocities,  each  of  which 

was  related  to  the  slope  of  the  stress  strain  curve  at  the  given  level  of  strain.  The  elastic 


stress  wave  would  travel  at  the  elastic  wave  speed  given  earlier  as, 


and  have  a 


magnitude  equal  to  O’  —  pcvo  .  The  plastic  wave  would  have  a  similarly  developed  wave 
velocity  of, 
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(18) 


where  Ep  is  the  slope  of  the  stress  strain  curve  in  the  plastic  region  and  cp  is  the  plastic 
wave  speed.  The  magnitude  of  the  stress  wave  would  then  be 


=  PCpvo 


(19) 


where  <jp  is  the  plastic  stress  wave  magnitude.  Figure  14  shows  the  stress  strain  curve  of  a 
bilinear  metal,  as  well  as  a  sample  wave  profile.  [6] 
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Figure  14  -  Bilinear  Stress-Strain  Curve  and  Corresponding  Wave  Profile 
An  alternative  analysis  did  not  use  the  bilinear  stress-strain  curve,  but  rather 

assumed  that  the  stress-strain  curve  was  concave  up  beyond  the  yield  stress,  see  Figure 

15.  The  plastic  wave  speed  was  just  shown  to  be 


c 


p 


(20) 


In  general  however,  when  the  slope  of  the  plastic  portion  of  the  stress-strain  curve  is  not 
linear,  the  plastic  wave  velocity  is  found  using 
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where  ^°/^£  is  the  slope  of  the  stress-strain  curve  at  a  given  strain.  This  plastic  velocity 

takes  the  form  of  the  particle  velocity  behind  the  shock.  A  material  with  a  concave-up 
stress-strain  curve  is  seen  in  Figure  15,  and  as  the  strain  is  increased  beyond  the  yield 
limit  at  point  A,  the  slope  of  the  stress-strain  curve,  and  therefore  the  velocity  of  the 
stress  wave  increases.  This  means  that  the  higher  stress,  increments  move  faster  than  the 
lower  stress  increments,  and  will  eventually  catch  up  to  the  lower  stress  increments,  at 
which  point  a  plastic  shock  front  is  formed,  as  shown  in  Figure  16. [6] 
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Figure  15  -  Concave-up  Stress-Strain  Curve 
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Figure  16  -  Shock  Formation 

The  rate-dependent  theory  was  developed  in  an  attempt  to  explain  some  of  the 
phenomena  that  could  not  be  described  by  the  rate  independent  theory.  The  first  attempt 
to  include  strain  rate  effects  was  called  the  overstress  model,  and  was  proposed  by 
Malvern,  with  the  form 


(7  = 


/(<?)+  In 


1  +  bs , 


V 


(22) 


where  f(s)  is  the  stress  found  from  a  quasi-static  stress-strain  curve,  h  is  a  constant,  and 


s  is  the  plastic  strain  rate.  An  alternative  form  of  the  overstress  model  provides  the 
plastic  strain  rate, 


(23) 
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where  the  term  cr  -  f{s )  is  the  overstress.  The  constants  a  and  b  are  material  specific 
constants,  that  are  used  to  correlate  with  experimental  data.  The  overstress  is  the 
difference  between  the  stress  applied  and  the  stress  that  would  occur  theoretically  in  a 
quasi-static  test  at  the  same  strain.  [6] 


2.2.3  Shock  Waves 

In  the  case  where  the  impact  velocity  is  much  greater  than  the  sound  speed  of  the 
material,  a  series  of  waves  led  by  an  elastic  wave,  which  is  overcome  by  a  plastic  wave, 
leads  to  the  formation  of  a  shock  wave.  Shock  waves  are  very  narrow  regions  in  a 
continuum  in  which  velocity,  temperature,  and  density  vary  in  a  nearly  discontinuous 
fashion.  This  discontinuity  causes  problems  with  the  conservation  equations  when  they 
are  in  a  differential  form,  because  all  properties  are  assumed  to  vary  in  a  continuous 
manner,  which  is  untrue  in  the  vicinity  of  a  shockwave. 

A  shock  wave  is  formed  by  the  coalescing  of  a  wave  front  of  various  speeds  and 
stress  levels  into  a  single  sharp  wave  front.  The  shock  wave  equations  are  developed  by 
a  simple  application  of  the  conservation  equations  presented  in  Section  2.1.  The  easiest 
case  to  study  is  the  case  where  a  shockwave  is  traveling  into  a  material  that  is  stationary, 
stress  free,  and  has  no  internal  energy  associated  with  it,  see  Figure  17.  For  simplicity,  it 
is  assumed  that  the  only  action  that  occurs  is  in  the  direction  of  the  shock’s  velocity.  This 
simplification  assumes  that  lateral  effects  can  be  ignored.  In  the  case  of  a  finite  element 
of  material,  this  assumption  may  not  be  accurate,  as  shock  reflections  are  likely  to  occur, 
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which  changes  the  state  of  stress  from  one-dimensional  to  three-dimensional.  This 
formulation  is  analogous  to  the  right  end  of  the  rod  in  Figure  5,  because  the  rod  is 
initially  stress  free  and  stationary.  The  two  states  identified  by  subscripts  0  and  1  in 
Figure  17  represent  the  physical  state  ahead  of  and  behind  the  shock,  which  is  shown 
traveling  at  a  velocity  of  U,  while  the  local  particle  velocity  is  given  by  u. 


The  conservation  of  mass  applied  across  the  shock  states  that  the  mass  flow 
entering  the  shock  must  equal  the  mass  flow  leaving  the  shock.  In  this  case,  the  frame  of 
reference  is  taken  to  be  moving  with  the  shock  at  a  velocity  of  U,  as  seen  in  Figure  18. 
The  velocity  of  the  particles  in  front  of  the  shock  was  assumed  to  be  zero,  which  is  why 
u0  is  set  equal  to  zero. 

Mathematically  the  conservation  of  mass  is  given  by, 

podA(uo  - U)At  +  p{dA(U  - ux)At  =  0  (24) 

where  the  mass  entering  the  shock  from  the  right  is  found  to  be  podA(ua  -U)At  since 
the  volume  of  material  moved  is  equal  to  the  area  through  which  the  shock  passes,  dA, 
multiplied  by  the  distance  the  material  moves  relative  to  the  shock,  which  is  (u0-U)At . 
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Similarly,  the  mass  entering  the  shock  from  the  left  is  found  to  be  pxdA{U  -  ul  )A t . 
Eliminating  the  differential  area  dA  and  time  At  from  both  sides,  and  remembering  that 
the  velocity  in  front  of  the  shock  is  zero,  we  obtain  the  first  of  the  shock  wave  equations, 

Pau  =  Pl(u  ~u i).  (25) 
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Figure  18  -  Moving  Coordinate  Frame 


The  conservation  of  momentum  is  developed  in  a  similar  manner  to  the  previous 
discussion  of  elastic  stress  waves,  except  that  the  new  terminology  is  used.  The  change 
in  momentum  across  the  shock  must  equal  the  impulse  applied,  or 


podAUAtux  =  a^dAAt 


(26) 


which  can  be  simplified  to  give  the  analogous  shock  wave  physics  definition  of  stress, 


(7,  =pJJux. 


(27) 
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The  last  of  the  conservation  equations  to  be  investigated  for  the  shock  wave  is  the 
conservation  of  energy.  The  conservation  of  energy  equation  says  that  initial  internal 
energy  plus  any  work  done  on  the  mass  is  equal  to  the  final  internal  energy, 

IE0  +  W0^  =  IEX  (28) 

where  IEq  is  the  initial  internal  energy,  IE]  is  the  final  internal  energy,  and  is  the 

work  to  go  from  state  zero  to  state  one.  The  internal  energy  is  a  combination  of  the 
internal  energy  source  per  unit  mass,  E,  plus  the  kinetic  energy  of  the  mass. 

IE0=E0+KE0  IE{  =  E{+  KEX .  (29) 

where  Eq  is  the  internal  energy  source  per  unit  mass  at  state  zero,  E/  is  the  internal  energy 
source  per  unit  mass  at  state  one,  KE0  is  the  initial  kinetic  energy,  and  KEi  is  the  initial 
kinetic  energy.  The  internal  energy  source  per  unit  mass  can  be  a  combination  of  things, 
such  as  a  chemical  reaction  that  releases  energy,  or  a  material  with  some  strain  energy 
that  is  stored.  Initially  the  kinetic  energy  is  zero,  because  u0  is  zero,  but  in  the  final  state 
the  kinetic  energy  is, 

KE,  =  ^(p(>UAt)(u,)2  (30) 

where  poUAt  is  the  mass  of  material  which  is  moving,  and  U]  is  the  velocity  at  which  it 
is  moving.  After  the  shock  has  passed,  the  internal  energy  source  is  found  to  be 

{Pi(U  -  ul)At)El  (31) 
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where  the  mass  is  found  using  the  velocity  relative  to  the  shock.  Combined,  the  previous 
two  equations  give  the  internal  energy  after  the  shock  has  passed  as, 

IE,  =  i  (pJJhtXu,  f  +  (A  (U  -  u,  )At)El  (32) 

The  internal  energy  source  in  front  of  the  shock  is, 

IEo  =  C PoUAt)Eo  (33) 

where  the  term  in  the  parentheses  is  the  mass  of  material  with  the  internal  energy  Ea. 
Because  the  mass  in  front  of  the  shock  was  assumed  stationary,  this  also  happens  to  be 
the  internal  energy  of  the  mass,  since  the  kinetic  energy  is  zero. 

Lastly,  work  is  known  to  be  a  force  carried  out  over  some  distance,  which  using 
the  nomenclature  presented  here  is  shown  to  be, 

W0^t  =  crt(utAt)  (34) 

there  is  no  area  in  this  fonnulation,  because  all  properties  are  assumed  to  be  applicable 
over  a  common  area,  and  it  is  therefore  eliminated  at  the  outset. 

Combining  the  previous  three  equations  gives  the  conservation  of  energy  equation 
for  the  case  of  a  moving  shock  wave, 

(poUAt)Eo  +  oxux  At  =  (px  ( U  -  ux )  A  t)Ex  +  i  (p0UAt)(ux  f  (35) 

eliminating  At  and  rearranging  so  that  the  internal  energy  terms  are  on  the  same  side, 
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(a  (.U  -  U>  ))E 1  -  PoUEo  =  A“l  -  ^  (A^XA  ): 


<-My 


The  second  term  on  the  left  side  of  the  equation  can  be  simplified  using  the  conservation 
of  mass,  and  the  second  term  on  the  right  side  can  be  simplified  using  the  conservation  of 
momentum,  giving 

(A  (U  ~  ui  ))Ei  ~Pi(U-  Ux  )E„  =  CT,M,  -  ^  <7,1 /, 

Simplifying  the  right  side  and  dividing  through  by  px(U  -ux)  gives 


E,~E .= 


At^-Mi)' 


Solving  the  conservation  of  momentum  equation  for  ui  and  substituting  into  the  above 
equation  gives 


EX~E0 


1 

~G\ 

2 


f 


V 


U(Pl-Po) 

A 


A 


y 


U(Pl~Po)^ 


A 


V 


A 


y 


Eliminating  £/  and  finding  a  common  denominator  on  the  bottom  gives 


^  CTi  (a  -  Po  ) 


A  A 


<4<» 
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Simplifying  one  step  further  gives  the  commonly  used  conservation  of  energy  equation 
for  shock  waves, 


E>-E„=\ct , 


f 


1  1 


\Po  Pi) 


(41) 


This  equation  shows  that  the  change  in  energy  between  the  two  states  is  a  function  of  the 
applied  stress  multiplied  by  the  difference  between  the  reciprocals  of  the  density.  The 
conservation  equations  by  themselves  present  an  incomplete  picture  of  the  physical  state 
around  the  shock  wave.  Investigation  of  the  three  conservation  equations  shows  that  they 
contain  a  total  of  five  unknowns,  see  Table  1. 


Law 

Table  1  -  Conservation  Equations 
Equation 

Unknowns 

Conservation 

PoU  =  PliU-u,) 

p,  U,  u 

of  Mass 

Conservation 

=PoUui 

p,  U,  U,  (7 

of  Momentum 

Conservation 
of  Energy 


Et-E0=la, 


f 


\ _ 1_ 

kPo  Px) 


p,  (7,  E 


Unknowns  U,  u,  E,  a,  p 


Having  five  unknowns  but  only  three  equations  to  describe  their  interaction 
presents  an  obvious  problem,  which  is  solved  through  the  use  of  a  Hugoniot  curve  or  an 
equation  of  state.  A  Hugoniot  curve  is  the  locus  of  all  attainable  shock  states  that  are 
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possible  in  a  material,  and  is  often  presented  in  a  form  relating  pressure  and  volume.  The 
Hugoniot  curve  is  similar  to  a  stress-strain  curve  in  uniaxial  stress,  except  it  relates 
pressure  and  volume  in  a  hydrostatic  loading  situation. 

A  Hugoniot  is  different  from  a  stress-strain  curve  in  that  it  is  not  developed  from 
one  experiment  that  follows  the  loading  path,  as  a  stress-strain  curve  does.  Instead,  a 
Hugoniot  curve  is  developed  using  a  large  number  of  planar  impact  experiments  to 
describe  the  relationship  between  the  hydrostatic  pressure  and  specific  volume.  Since 
each  point  in  a  Hugoniot  curve  represents  a  separate  experiment,  a  Hugoniot  curve  does 
not  represent  a  loading  path  that  is  followed,  but  rather  each  point  is  an  equilibrium  point 
for  a  specific  experiment.  A  generic  Hugoniot  curve  is  seen  in  Figure  19. [7] 


Figure  19  -  Hugoniot  Curve  Showing  Loading  Path  and  Unloading  Path 
On  this  plot,  the  Hugoniot  curve  is  marked  with  an  H.  When  an  impact  occurs 

with  an  initial  velocity  of  va,  the  loading  path  follows  the  line  from  point  A  where  the 

material  starts  out  with  zero  pressure  but  a  high  velocity,  to  point  B  along  what  is  called 
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the  Rayleigh  line.  Thus,  it  is  seen  that  the  loading  does  not  follow  the  Hugoniot,  but 
rather  occurs  along  a  straight  line  connecting  the  initial  state  with  the  point  along  the 
Hugoniot  curve  relating  to  the  peak  pressure  of  impact,  which  is  marked  here  as  Ph- 
Unloading  in  this  case  occurs  isentropically  along  the  line  marked  S,  which  is  also  not 
along  the  Hugoniot  curve. [7] 

Since  the  Hugoniot  curves  are  developed  under  uniaxial  strain  shock  wave 
conditions,  they  are  only  valid  in  certain  restricted  situations,  which  leads  to  the  use  of 
equations  of  state,  which  are  more  general. [5]  In  most  computer  codes  that  solve  impact 
problems,  an  equation  of  state  is  used  to  relate  internal  energy,  pressure,  and  volume.  An 
equation  of  state  is  developed  using  planar  impact  experiments,  which  can  be  used  to 
develop  the  Hugoniot,  and  curve-fitting.  Examples  of  equations  of  state  will  be  presented 
in  a  later  section. 

As  was  mentioned,  a  shock  wave  is  formed  when  higher  velocity  wavelets 
overtake  slower  velocity  wavelets,  which  leads  to  an  instantaneous  change  in  the  state  of 
the  material.  Because  a  shock  wave  requires  waves  of  different  velocity  to  exist  in  a 
material,  it  is  impossible  for  an  elastic  shock  wave  to  be  formed,  because  the  elastic 
sound  speed  of  a  material  is  constant. 

A  shock  wave  might  look  similar  to  an  elastic  wave  in  that  they  each  have  a  very 
definite  wave  front,  but  the  difference  is  that  in  an  elastic  wave,  the  wave  front  is  due  to 
the  impact,  and  never  changes,  while  in  a  shock  wave,  the  wave  front  develops  over  a 
short  period  of  time,  and  can  have  a  varying  amplitude. [5] 
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2.3  Constitutive  Equation 


The  relationship  between  stress  and  strain  in  continuum  mechanics  codes  is 
dictated  by  a  constitutive  equation.  In  most  finite  element  codes,  stress  is  assumed  to  be 
quasi-static,  which  means  that  the  loading  is  applied  so  slowly  that  there  aren’t  any 
dynamic  loading  effects.  In  quasi-static  cases,  the  most  common  constitutive  equation 
used  is  the  classic  Hooke’s  Law  equation, 

cr  =  f(s,Em)  (42) 


where  a  is  the  stress,  s  is  the  strain,  and  Em  is  the  modulus  of  elasticity.  In  many 
situations  however,  it  is  inappropriate  to  assume  that  stress  is  applied  quasi-statically, 
because  of  this,  Hooke’s  Law  will  only  be  used  in  cases  where  the  stress  is  below  the 
yield  stress  of  the  material.  In  cases  where  the  applied  stress  is  greater  than  the  yield 
stress,  it  is  necessary  to  account  for  dynamic  loading  effects.  The  most  common  way  to 
account  for  dynamics  in  a  continuum  mechanics  problem  is  to  include  strain  rate  as  a 
variable  in  the  constitutive  equation.  In  general  this  becomes, 


f 


°  =  f 


£>  £’  Em 


V 


(43) 


where  a  is  the  strain  rate  applied.  In  some  cases,  constitutive  equations  will  also  be  a 
function  of  internal  energy  and  damage.  [8] 
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CTH  provides  numerous  constitutive  equations  with  which  stress-strain  behavior 
can  be  modeled.  Most  of  these  equations  will  be  of  little  use  in  this  impact  study, 
because  constitutive  equations  tend  to  be  very  problem  specific.  Constitutive  equations 
exist  for  metals,  ceramics,  concrete,  and  soil  amongst  others.  For  this  problem,  only 
those  equations  dealing  with  metals  are  applicable. 

One  of  the  most  basic,  yet  still  valuable,  constitutive  models  available  in  CTH  is 
the  Johnson-Cook  Strength  Model.  This  model  presents  the  von  Mises  flow  stress  as 


cr  =  (A  +  Bsn  )(1  +  C  In  £■  )(1  -  T*m ) 


(44) 


where  a  is  the  von  Mises  flow  stress,  s  is  the  equivalent  plastic  strain,  s  P  is  the  plastic 

strain  rate  nonnalized  by  a  strain  rate  of  1.0s'1,  T*  is  defined  below,  and  A,  B,  C,  m,  and 
n  are  the  Johnson-Cook  coefficients  for  the  given  material.  The  Johnson-Cook 
viscoplastic  material  model  accounts  for  temperature  via  the  homologous  temperature,  T* 
which  is  given  as 


rj-i  =1 


T-T. 


T  -T 

melt  room 


(45) 


where  T  is  the  absolute  temperature,  Troom  is  the  ambient  temperature,  and  Tmeit,  is  the 
melting  temperature  of  the  material.  [5] 

There  are  two  minor  disadvantages  to  the  Johnson-Cook  model.  The  first  is  that  it 
presents  strain  rate  sensitivity  as  being  independent  of  temperature,  which  in  general  is 
not  the  case.  However,  by  keeping  strain,  strain  rate,  and  temperature  uncoupled,  it 
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becomes  relatively  straightforward  to  determine  the  Johnson-Cook  coefficients  from  a 
few  simple  experiments  at  various  temperatures  and  strain  rates.  The  second 
disadvantage  of  the  Johnson-Cook  model  is  that  it  is  strictly  a  mathematical  curve-fit  of 
experimental  data,  and  is  therefore  not  built  upon  a  base  of  physics. 

One  constitutive  model  that  overcomes  the  disadvantages  of  the  Johnson-Cook 
model  is  the  Zerilli- Armstrong  model.  This  model  accounts  for  the  interdependence  of 
strain,  strain  rate,  and  temperature,  and  does  it  from  a  dislocation  dynamics  basis. 
Because  it  is  based  on  dislocation  dynamics,  it  takes  different  forms  for  different  metal 
structures.  The  face-centered-cubic  form  of  the  Zerilli-Annstrong  model  is 


C7  =  c0+c,4e  ■e^T+c‘TAnc 


(46) 


where  the  C,’s  are  the  Zerilli-Annstrong  coefficients,  s  is  the  equivalent  plastic  strain,  s 
is  the  equivalent  strain  rate,  and  T  is  the  absolute  temperature.  The  body-centered-cubic 
form  is 


<r  =  C0  +  C;  •  e-cir+c*r  'ne  +  cse' 


(47) 


where  n  is  a  strain-hardening  coefficient.  Simulations  using  the  Zerilli-Annstrong  and 
Johnson-Cook  models  show  slightly  better  results  for  the  Zerilli-Annstrong  model  when 
compared  to  experimental  results;  however  neither  model  was  especially  effective  for 
very  large  strains.  [9] 
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A  third  constitutive  model,  which  also  utilizes  dislocation  dynamics,  is  the 
Steinberg-Guinan-Lund  Model.  This  model  defines  the  dynamic  yield  stress,  the  shear 


modulus,  and  also  the  melting  temperature.  The  dynamic  yield  stress  is  given  by 


g(p,t) 


G 


(48) 


where  T?  is  the  thermally  activated  component,  Yj  is  the  yield  stress  at  the  Hugoniot 

elastic  limit,  sp  is  the  plastic  strain,  sp  is  the  plastic  strain  rate,  G0  is  the  initial  shear 
modulus,  and  G(P,T)  is  the  shear  modulus,  found  using 


g(p,t)=g0 


,  AP 
+  G 


B(T-  300) 


(49) 


where  A  and  B  are  material  constants,  P  is  the  pressure,  p  is  the  density  ratio  of 


and  T  is  the  temperature  in  Kelvin.  The  function  f(ep )  is  the  work  hardening  function, 
which  is  found  using 


ASp)=^+P(£P  +£i)\  GO) 


where  /?  and  S\  are  material  constants.  Finally,  the  thermally  activated  component  is  a 
result  of  the  plastic  strain  rate, 
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(51) 
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where  Yp  is  the  Peierl’s  stress,  2Uk  is  the  energy  required  to  form  two  kinks  in  a 
dislocation  segment,  and  Cj  and  C?  are  material  constants. 

Melting  of  material  is  modeled  in  the  Steinberg-Guinan-Lund  model  as  well.  In 
this  model,  the  melting  temperature  for  a  given  state  is  found  via 


Mil 

2  a\ 

l  n). 

Aro-a-%) 


m,o 


(52) 


where  Tm>0  is  the  melting  temperature  at  constant  volume  and  a  and  y0  are  material 
constants.  A  table  of  the  constants  used  in  the  specific  problem  studied  are  presented  in 
Table  2.  The  materials  that  were  actually  used  in  CTH  are  iron  for  1080  steel  and 
VascoMax  250  for  VascoMax  300.  The  yield  stress  and  density  for  iron  and  VascoMax 
250  were  changed  to  match  the  values  of  1080  steel  and  VascoMax  300. [2] 
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Table  2  -  Material  Constants  used  in  Constitutive  Models 


Johnson-Cook  Constants 

1080  Steel  (modeled  as  iron) 

Steinberg-Guinan-Lund  Constants 

VascoMax  300  (modeled  as  VascoMax  250) 

Constant 

Value 

Constant 

Value 

A 

1.7526  x  109 

A 

2.06  x  10'12 

B 

3.8019  x  109 

B 

3.15  x  10'4 

C 

0.06 

Po 

8.129 

m 

0.55 

G0 

7.18  x  1011 

n 

0.32 

Y0 

1.447  x  1010 

Tm 

1835.7 

y 

1  max 

2.5  x  1010 

Tm,0 

2310 

n 

0.5 

a 

1.2 

P 

2.0 

Yo 

1.67 

2.4  Equation  of  State 

It  is  common  when  solving  dynamic  mechanics  problems  to  break  down  stress 
and  strain  into  two  components,  the  hydrostatic  or  volumetric  stress  or  strain  and  the 
deviatoric  stress  or  strain, 

H=hJ+ki]  (S3) 

where  [a]  is  the  stress  tensor,  [ah]  is  the  hydrostatic  stress  tensor,  and  [a<i]  is  the 
deviatoric  stress  tensor.  The  hydrostatic  stress  is  often  called  the  volumetric  stress 
because  it  is  the  stress  that  develops  a  volume  change  for  a  given  parallelepiped  of 
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material,  while  the  deviatoric  stress  is  associated  with  a  shape  change.  In  impact 
problems,  these  two  varieties  of  stress  are  handled  via  two  separate  relationships.  The 
first  of  these  relationships,  the  deviatoric  stress,  was  discussed  in  the  previous  section  on 
constitutive  equations.  The  second  relationship  deals  with  the  hydrostatic  stress,  and  is 
called  the  equation  of  state.  The  two  are  taken  separately  because  it  has  been  found  that 
hydrostatic  stress  is  virtually  independent  of  strength  and  plasticity,  while  deviatoric 
stress  is  only  slightly  dependent  upon  pressure.  [5]  Additionally,  equations  of  state  are 
needed  to  model  how  pressure,  density,  and  energy  relate  when  compressibility  effects, 
and  irreversible  processes  such  as  shock  waves  are  included  in  the  problem.  [8] 

The  equation  of  state  of  a  material  describes  the  relationship  between  pressure, 
specific  volume,  and  internal  energy,  and  can  be  shown  in  a  general  form  by 

E  =  E(P,  V)  <}4) 

where  E  is  the  internal  energy,  P  is  the  pressure,  and  V  is  the  specific  volume.  An 
alternative  form,  often  used  in  computer  codes  is, 

P  =  P(p,E)  (55) 

The  Mie-Griineisen  equation  of  state  is  a  simple  equation  of  state  that  is  very 
good  for  modeling  high-pressure  shock  related  events. [10]  The  Mie-Griineisen  equation 
of  state  is  based  upon  statistical  mechanics,  using  the  energy  of  individual  atoms  to  arrive 
at  thermodynamic  equations.  The  Hugoniot  pressure  is  used  as  a  baseline  in  the  Mie- 
Griineisen  equation  of  state  and  is  given  by, 
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PH  —  Cj  n  +  C0  /u  +  C^ju 


(56) 


where  Ph  is  the  Hugoniot  pressure,  the  C,  ’v  are  constants,  and  //  is 


P_ 

Po 


-1 


f57; 


The  C  parameters  in  the  equation  for  the  Hugoniot  pressure  are  only  for  a  case  where 
density  increases.  If  density  decreases,  C2  and  C?  are  zero.  The  pressure  is  then 
calculated  with 


P  =  P 


f 


H 


1- 


Y  ju 


+  Tp{E-E„) 


(58) 


where  E  is  the  internal  energy  per  unit  mass,  E0  is  the  internal  energy  per  unit  mass  at 
ambient  conditions,  and  T  is  a  constant  called  the  Griineisen  parameter.  The  Griineisen 
parameter  is  assumed  to  be  independent  of  temperature  and  only  a  function  of  specific 
volume,  and  is  represented  as 


r  =  v 


KdEj 


(59) 


Another  equation  of  state  often  used  is  the  Tillotson  equation  of  state.  This 
equation  also  uses  a  quadratic  approximation  for  pressure,  and  also  has  different  versions 
for  different  density  values.  When  the  density  is  greater  than  the  ambient  density  and  the 
internal  energy  is  less  than  the  sublimation  energy,  the  pressure  is  given  by 
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(60) 


P  —  Pn  +  A/u  +  Bp 


where  A  and  B  are  material  constants,  and  P,T  is, 


P*=Ep 


a  +  - 


+  1 


(61) 


where  a  and  b  are  constants,  and  77  =  — . 

Po 

When  density  is  less  than  ambient  and  internal  energy  is  greater  than  sublimation 
energy,  the  pressure  is  given  by, 


P  -  ciEp 


bEp 


+  Ape 


-P 


>1- 1 


1  + 


EoW 


;;-l 


(62) 


where  a  and  /?  are  constants.  [2]  These  are  just  two  of  the  many  different  equations  of 
state  that  are  commonly  used  in  solving  impact  problems,  others  exist  for  gases  and 
explosives,  and  are  available  for  use  in  CTH  if  needed. 

The  equation  of  state  used  in  this  investigation  isn’t  actually  an  equation  at  all.  It 
is  in  fact  simply  a  table  that  correlates  pressure,  energy,  and  density  at  various  states.  In 
CTH,  this  equation  of  state  is  called  the  SESAME  model.  Two  major  advantages  of  a 
tabular  equation  of  state  are  that  there  is  no  need  to  calculate  equation  of  state  variables, 
as  they  are  simply  part  of  a  table,  and  that  by  using  a  tabular  equation  of  state  the  exact 
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physical  state  is  used  as  opposed  to  an  assumed  state,  i.e.  a  quadratic  form  as  in  the  Mie- 
Griineisen  and  Tillotson  equations  of  state.  This  can  be  very  important  if  the  pressures 
applied  are  high  enough  that  a  material  will  change  state  from  solid  to  liquid  or  liquid  to 
gas. 


2.5  Artificial  Viscosity 

In  a  continuum,  it  is  assumed  that  physical  properties  within  the  continuum  vary 
in  a  smooth  manner,  meaning  that  there  are  no  instantaneous  jumps  in  properties.  This  is 
a  requirement  of  continuum  mechanics  codes,  which  solve  the  conservation  equations  in 
differential  form.  The  problem  that  is  encountered  in  high  velocity  impact  problems  is 
that  at  high  velocities,  a  shock  wave  is  formed.  Shock  waves  are  in  effect  instantaneous 
jumps  in  pressure,  density,  and  velocity,  which  make  it  impossible  to  solve  the 
conservation  equations  across  a  finite  mesh  size.  This  means  that  if  a  shock  wave  occurs 
within  an  element  of  a  mesh,  it  may  not  even  be  observed.  The  most  common  method  of 
eliminating  this  problem  is  to  spread  the  shock  wave  across  a  few  elements,  eliminating 
the  discontinuity  that  causes  problems  for  the  differential  equations. 

The  method  of  spreading  shocks  across  a  number  of  cells  used  in  most  continuum 
mechanics  codes,  including  CTH,  is  artificial  viscosity.  Artificial  viscosity  was 
developed  based  upon  viscosity  used  in  fluid  mechanics.  The  one-dimensional 
conservation  of  momentum  equation  for  fluid  dynamics  contains  a  viscous  term, 
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d  (  dv ^ 

u  — 

V  dx) 


dx 


(63) 


where  //  is  the  viscosity  of  the  fluid.  This  led  to  the  idea  of  adding  a  similar  term  to  the 
pressure  tenn  in  the  conservation  equations 


c '{p  +  q ) 

dx 


(64) 


where  q  is  the  artificial  viscosity,  which  is  of  the  same  form  as  the  viscosity  used  in  fluid 
mechanics.  It  was  found  that  the  best  way  to  smear  the  shock,  without  affecting  the 
solution  away  from  the  shock  was  to  add  a  linear  term  and  a  quadratic  term.  These  two  q 
tenns  were  found  to  be, 
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q  2  =  (b2Ax)cp  — 
dx 


(65) 


where  bi  and  />2  are  constants  used  to  control  how  the  shock  is  spread  across  zones.  An 
example  of  the  use  of  artificial  viscosity  is  seen  in  Figure  20.  The  first  graph  shows  a 
theoretical  shock,  dashed  line,  and  the  simulated  behavior  that  would  be  seen  without 
applying  artificial  viscosity  to  the  problem.  The  second  graph  shows  the  same  shock, 
except  in  this  case,  artificial  viscosity  is  included.  It  is  obvious  that  artificial  viscosity 
greatly  decreases  the  oscillations  that  would  be  simulated.  By  increasing  the  fineness  of 
the  mesh,  it  would  be  possible  to  improve  the  shock  front  even  more,  to  more  closely 
simulate  a  true  shock  wave.  [8] 
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Chapter  3  -  Methodology 


3.1  CTH  Solution  Method 

CTH  is  a  software  system  developed  by  Sandia  National  Laboratories  in  the  late 
1980’s  and  early  1990’s,  designed  to  “model  multidimensional,  multi-material,  large 
deformation,  strong  shock  wave  physics”  [11]  problems.  CTH  is  a  hydrocode  computer 
program,  which  evolved  from  hydrodynamic  codes,  which  treated  materials  as  being 
fluid  like,  and  therefore  possessing  no  strength  in  tension  or  shear.  CTH  however  is  an 
advanced  hydrocode,  which  does  have  the  capability  of  modeling  strength  in  materials. 

CTH  differs  from  other  continuum  mechanics  codes  in  the  solution  scheme  that  it 
uses.  In  general  there  are  two  different  ways  that  a  continuum  can  be  described,  the 
Lagrangian  or  material  description  and  the  Eulerian  or  spatial  description.  The 
Lagrangian  or  material  description  of  a  continuum  essentially  takes  the  material  to  be 
acted  upon  and  divides  it  up  into  smaller  pieces  and  solves  the  conservation  equations  by 
following  specific  pieces  of  material.  The  Eulerian  or  spatial  description  defines  a 
volume  in  space  and  solves  the  conservation  equations  by  tracking  what  goes  through  this 
volume. 

There  are  advantages  and  disadvantages  for  using  either  of  these  two  meshing 
alternatives.  In  a  Lagrangian  mesh,  it  is  very  easy  to  define  boundaries  and  to  calculate 
properties  for  a  specific  piece  of  material,  but  it  is  not  a  very  good  choice  if  one  wishes  to 
solve  a  problem  where  very  large  defonnations  occur.  Large  deformations  in  a 
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Lagrangian  mesh  tend  to  deform  the  elements  far  too  much  to  provide  accurate  answers. 
In  an  Eulerian  mesh,  large  deformations  are  not  a  problem,  since  material  is  not  tracked 
specifically,  and  therefore  it  can  deform  as  much  as  needed,  without  causing  numerical 
problems.  Eulerian  meshes  however,  tend  to  make  it  more  difficult  to  define  material 
boundaries  when  mixing  occurs,  and  provide  no  ready  way  to  watch  what  happens  to  a 
specific  piece  of  material. 

In  problems  involving  high  velocity  impact,  it  is  seen  that  neither  a  Lagrangian 
solution  nor  an  Eulerian  solution  is  ideal,  because  in  general  large  deformations  are 
exactly  what  is  to  be  studied,  yet  it  is  still  desired  to  see  physical  properties  related  to  a 
specific  piece  of  the  material.  Ideally,  the  best  properties  of  both  Lagrangian  and 
Eulerian  solutions  would  be  combined  into  one  package  to  solve  high  velocity  impact 
problems.  This  is  exactly  what  CTH  has  been  designed  to  do. 

The  method  CTH  uses  to  solve  impact  problems  is  known  as  a  two-step  Eulerian 
solution  scheme.  The  CTH  solution  starts  with  an  Eulerian  mesh,  which  during  the  first 
step  is  allowed  to  defonn  in  a  Lagrangian  manner.  The  second  step  takes  the  deformed 
mesh  and  maps  it  back  to  the  original  Eulerian  mesh.  CTH  provides  for  one-dimensional 
solutions,  two-dimensional  rectangular  and  cylindrical  solutions,  and  three-dimensional 
rectangular  solutions. 

The  quantities  to  be  calculated  in  a  CTH  solution,  except  for  velocity,  are  all 
assumed  to  be  constant  across  each  individual  cell,  and  centered  within  the  cell.  The 
velocity  of  material  is  assumed  to  act  on  the  cell  face.  In  order  to  solve  the  conservation 
of  momentum  equation,  it  is  necessary  to  create  a  staggered  cell  structure,  where  the 
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edges  of  the  staggered  cell  are  located  at  the  center  of  the  original  cell.  Since  all 
quantities  except  for  velocity  are  centered  over  the  original  mesh,  these  quantities  are  all 
considered  to  be  based  in  time  on  the  whole  time  step,  while  velocity  is  based  on  a  half 
time  step,  because  it  is  spatially  on  a  staggered  cell.  An  example  of  the  location  of  the 
cell  quantities  and  a  depiction  of  the  staggered  cell  structure  are  seen  in  Figure  21.  [11] 
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Figure  21  -  CTH  Cell  Structure 

As  was  mentioned  earlier,  CTH  utilizes  a  two-step  solution  scheme.  The  first  step 
integrates  the  conservation  equations,  which  are  presented  as  a  finite  volume 
approximation  across  the  time  step.  The  integration  of  the  conservation  equations  is  done 
explicitly,  which  is  simpler  than  implicitly,  but  requires  a  smaller  time  step,  to  ensure 
stability.  The  user  does  not  control  the  time  step  that  CTH  utilizes  directly,  since  it  is 
calculated  to  ensure  that  a  wave  will  not  skip  a  cell  during  one  time  step.  Since  this  step 
is  a  Lagrangian  step,  with  the  mesh  attached  to  specific  pieces  of  material,  it  is  easy  to 
see  that  there  is  no  mass  that  moves  out  of  the  mesh,  and  therefore  for  the  first  step,  the 
conservation  of  mass  equation  is  satisfied  automatically.  [11] 

The  remap  step  utilized  by  CTH  returns  the  defonned  mesh  to  its  original  state, 
and  calculates  the  amount  of  volume,  mass,  momentum,  and  energy  that  must  be  moved 
into  each  cell.  Volume  flux  between  neighboring  cells  is  calculated,  and  a  high- 


original  Mesh 
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resolution  interface  tracking  subroutine  is  used  to  determine  which  material  moves.  The 
next  step  calculates  how  much  mass  and  internal  energy  moves  into  the  cells.  Lastly, 
momentum  and  kinetic  energy  are  moved,  because  these  two  are  dependent  upon  the 
mass.  [11] 

In  the  vast  majority  of  the  structure,  the  remap  step  will  likely  be  rather  simple, 
because  only  one  material  is  in  a  given  area,  and  therefore  it  isn’t  necessary  to  detennine 
how  much  of  which  material  moves  where.  However  in  areas  where  there  are  multiple 
materials,  the  remap  step  takes  on  added  difficulties,  because  it  is  necessary  to  keep  track 
of  multiple  individual  bits  of  material,  which  in  fact  might  be  only  a  very  small  amount 
of  material.  [11] 

In  the  case  where  multiple  materials  are  present  in  a  cell,  the  two-dimensional 
solution  approximates  the  material  boundaries  as  straight  lines.  The  volume  fraction  of 
each  corner  is  found  by  averaging  the  volume  fractions  of  the  four  cells  surrounding  the 
comer  point.  The  values  between  these  corner  values  are  assumed  to  vary  linearly,  and 
can  be  interpolated  easily.  The  equivalent  volume  fractions  on  adjacent  sides  are 
connected  via  a  straight  line,  see  Figure  22.  For  each  case  of  corresponding  volume 
fractions,  i.e.  0.5  and  0.5,  the  volume  below  the  diagonal  line  is  calculated.  This 
calculated  value  is  compared  to  the  cell  volume  fraction,  and  if  they  are  not  equal,  a 
different  pair  of  volume  fractions  is  used,  until  the  calculated  value  and  the  cell  value  are 
equal.  In  Figure  22,  the  volume  flux  calculated  earlier  is  drawn  as  the  volume  to  the  right 
of  the  dashed  line.  Using  these  two  lines,  the  material  to  be  fluxed  to  the  right  is  the 
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volume  of  material  to  the  right  of  the  dashed  line  and  below  the  diagonal  line.  This  is 
repeated  for  all  four  sides  and  for  each  material  present  in  a  cell.  [11] 
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Figure  22  -  Material  Flux  Between  Cells 
A  difficulty  arises  when  CTH  attempts  to  conserve  mass,  momentum,  and  kinetic 

energy  during  the  remap  step.  The  problem  is  that  CTH  is  trying  to  conserve  three 

quantities,  while  only  using  two  variables  to  do  it.  CTH  provides  three  different  methods 

to  handle  this  problem.  The  first  method  chooses  to  conserve  momentum,  but  to  change 

kinetic  energy  into  internal  energy.  As  can  be  imagined,  this  case  could  cause  unrealistic 

temperature  increases  in  the  cell,  if  the  kinetic  energy  is  high.  The  second  method  is  the 

reverse  of  the  first,  conserving  kinetic  energy,  but  not  momentum.  This  method  does  not 

give  satisfactory  answers.  The  final  method  for  conserving  mass,  momentum,  and  kinetic 

energy  is  the  same  as  the  first  method,  except  in  the  case  when  the  receiving  cell  has 

momentum  of  a  different  sign  than  the  donating  cell.  In  this  special  case,  the  kinetic 

energy  is  added  to  the  internal  energy  of  the  receiving  cell.  [12] 

Since  multiple  materials  are  capable  of  occupying  a  given  cell,  it  is  necessary  to 

understand  how  pressure  and  temperatures  are  handled  within  the  cell  for  each  of  these 

materials.  CTH  has  the  capability  of  having  multiple  temperatures  and  pressures  within  a 

given  cell.  The  first  option,  allows  each  material  to  have  separate  temperatures  and 
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pressures,  but  does  not  allow  pressure  relaxation,  which  causes  problems  with  the 
conservation  equations.  This  option  distributes  volume  and  energy  based  upon  the 
material  volume  fractions  within  each  cell,  and  can  cause  problems  when  materials  have 
compressibility  mismatches.  The  second  option  is  similar  to  the  first  in  its  treatment  of 
mixed  cells,  however  it  distributes  volume  and  energy  based  upon  the  volume  fraction 
cubed  divided  by  the  material  mass.  This  option  seems  to  provide  better  results  than  the 
first  option  in  most  cases.  The  third  option  allocates  work  done  in  cells  based  upon 
compressibility  of  the  given  materials,  and  allows  pressure  relaxation.  This  option  gives 
better  results  than  either  the  first  or  second  option,  but  the  algorithm  is  less  robust.  The 
last  option  allocates  work  done  based  upon  volume  fraction,  but  does  so  in  a  way  that 
prevents  materials  with  very  low  volume  fractions  from  changing  volume.  This  option 
also  allows  pressure  relaxation.  [12] 

One  last  issue  for  mixed  cells  needs  to  be  discussed,  that  is  how  to  calculate  the 
yield  strength  in  a  cell  of  mixed  materials.  The  first  option  uses  the  volume  fraction 
averaged  yield  strength.  The  second  option  is  the  volume  averaged  yield  strength 
nonnalized  by  the  sum  of  the  volume  fractions  of  materials  that  support  a  shear  loading. 
The  last  option  is  that  the  strength  in  a  cell  of  mixed  materials  is  zero  unless  the  cell 
contains  a  material  and  a  void;  in  which  case,  a  volume  averaged  yield  strength  is 
used.  [12] 

One  defining  characteristic  of  CTH  when  compared  to  a  typical  hydrocode  is  that 
CTH  models  the  ability  of  materials  to  support  a  deviatoric,  or  shearing,  stress.  In 
hypervelocity  gouging  problems,  one  of  the  main  deformation  causing  processes  is 
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believed  to  be  the  sliding  of  the  slipper  along  the  rail.  CTH  provides  two  options  to 
model  a  sliding  interface  between  two  materials. 

The  first  option  was  developed  to  model  penetration  events,  and  is  characterized 
by  the  movement  of  the  sliding  process  into  what  is  deemed  the  “soft”  material,  away 
from  the  interface  of  the  materials,  which  improves  numerical  results.  Moving  the 
sliding  interface,  allows  the  development  of  deviatoric  stress  at  the  interface.  This 
algorithm  unfortunately  is  not  supported  for  all  geometries,  and  is  unable  to  be  utilized  in 
a  parallel  computing  environment,  which  therefore  makes  it  inconvenient  to  use  in  many 
situations  that  aren’t  axisymmetric,  or  that  are  very  computationally  intensive,  since  it 
must  be  used  on  a  single  processor.  [12] 

The  alternative  algorithm  takes  a  different  approach  to  handling  a  sliding  interface 
than  the  other  algorithm  does.  Instead  of  moving  the  sliding  interface  away  from  the 
material  interface  as  the  previous  algorithm  does,  this  algorithm  sets  the  deviatoric  stress 
at  the  material  boundary  to  zero.  This  option  effectively  turns  the  projectile’s  surface 
into  a  liquid,  since  it  is  unable  to  support  any  frictional  forces.  [12] 
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3.2  Buckingham  Pi  Theorem 


In  his  work  on  the  gouging  problem  for  the  HHSTT,  Szmerekovsky  utilized  the 
Buckingham  Pi  theorem  to  scale  his  model  appropriately  so  that  he  represented  the  actual 
physical  problem  as  accurately  as  possible,  but  on  a  smaller  scale.  This  scaling  was 
necessary  to  keep  the  simulation  computationally  efficient. 

The  Buckingham  Pi  theorem  is  a  method  of  dimensionally  analyzing  a  problem, 
that  allows  models  to  be  scaled,  yet  still  be  consistent  with  some  physical  problem.  The 
Buckingham  Pi  Theorem  is  a  method,  by  which  a  physical  problem  can  be  broken  down 
into  a  set  of  invariant  products  that  must  be  satisfied  for  a  system  to  remain  consistent. 
The  Buckingham  Pi  Theorem  states  that  if  a  physical  system  consists  of  a  number  of 
dimensioned  quantities  {qi}m  that  are  each  products  of  a  set  of  j  independent  fundamental 
dimensions  Lj,  then  the  physical  law  can  be  described  by 

f(ql,q2,q3,.".,qm)  =  0  (66) 


where  m  is  the  number  of  dimensioned  quantities  which  are  to  be  used  in  the 
Buckingham  Pi  analysis. [13] 

A  fundamental  dimension  is  a  quantity  that  is  used  to  describe  a  dimensioned 
quantity.  Szmerekovsky  utilized  mass,  length,  time,  and  temperature  as  the  four 
fundamental  dimensions  in  his  analysis  and  represented  them  by  M,  L,  T,  and  6 
respectively.  Using  these  four  fundamental  dimensions,  it  is  possible  to  describe  most 


physical  quantities,  such  as  velocity,  which  is  j  or  force  ^  .  It  is  not  necessary  to  use 
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these  four  fundamental  dimensions;  it  is  only  necessary  to  ensure  that  the  fundamental 
dimensions  alone  can  describe  all  dimensioned  quantities.  For  example  one  could  use 

velocity  ^  and  time  T  if  one  wishes.  This  might  be  useful  in  a  case  where  a  length 

might  not  be  easily  measurable,  but  the  velocity  and  time  can  be  observed  easily.  [2] 

As  was  just  mentioned  for  the  case  of  velocity  and  force,  it  is  possible  to  represent 
any  dimensioned  quantity  as  a  product  of  fundamental  quantities  raised  to  some  power, 


Vi  = 


J  d\  J  d2 

^1  ^2 


(67) 


where  q,  is  the  dimensioned  quantity,  Lj  is  the  fundamental  dimension,  and  dk  is  the 
power  the  fundamental  dimension  is  raised  to.  In  the  case  of  force  this  leads  to, 


LX=M 

l2=l 

L,=T 
dt  =1 
f/2  —  1 
cl,  =  -2 


(68) 


This  process  is  repeated  for  all  of  the  dimensioned  quantities  that  are  to  be  scaled. 
These  dimensioned  quantities  are  then  combined  to  form  the  invariant  Pi  quantities,  in  the 
following  manner 

n  =  (<?,)“'(<?, (69> 
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where  the  at  ’s  are  an  exponent  that  is  to  be  detennined.  Utilizing  the  equation  describing 


the  dimensioned  quantities  in  terms  of  the  fundamental  dimensions  gives, 


n  =  (z/iz/2---z/”)rii(z/iz2 


This  equation  is  rearranged  so  that  all  of  the  Li  quantities  are  together  and  so  on.  This 
can  be  simplified  by  combining  the  powers  of  the  fundamental  dimensions,  giving 


where  the  exponents  ft  can  be  described  as 


A 

dn 

dn 

Pi 

-  = 

d2l 

d22 

A 

_dn\ 

dn2 

dim 
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2m 


d. 


ax 

a2 


(71) 


(72) 


Since  this  system  must  be  applicable  for  any  physical  law,  a  function  of  the  Pi 
invariants  must  be  invariant.  This  in  turn  requires  that  the  exponents  /?bc  equal  to  zero, 
otherwise  when  a  fundamental  dimension  changed,  the  Pi  invariant  would  change.  By 
solving  the  above  equation,  it  is  possible  to  determine  which  products  of  dimensioned 
quantities  must  be  maintained  invariant.  An  example  follows  which  will  show  how  to 
perform  a  Buckingham  Pi  analysis  on  the  problem  of  a  cylindrical  projectile  impacting  a 
flat  plate  at  an  angle,  see  Figure  23. 
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Figure  23  -  Rod  Impact  Example 

The  dimensioned  quantities  to  be  utilized  follow  from  Szmerekovsky’s  work,  and 
are  shown  in  Table  3.  Additionally,  the  description  of  these  quantities  in  their 
fundamental  dimensions  is  provided. 

_ Table  3  -  Buckingham  Pi  Example _ 


Dimensioned  Quantity 

Symbol 

Fundamental  Dimensions 

density 

P 

ML'3 

diameter 

d 

L 

length 

1 

L 

velocity 

V 

LT'1 

material  speed  of  sound 

c 

LT'1 

internal  energy 

E 

ML2T'2 

energy  source 

S 

ML2T'2 

yield  strength  in  compression 

<Jy,c 

ML-1T'2 

elastic  modulus 

f 

‘-'in 

ML_1T'2 

shear  modulus 

G0 

ML_1T'2 

time 

t 

T 
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These  quantities  are  then  raised  to  a  power  a  and  multiplied  together  to  fonn  the 
invariant  II. 


n  =  (pf  (df‘  (/)“•  (v)“*  (cf’  (£)“«  (5)“'  (<x„  Y’  {. Em  )“>  (G„  ( t ) 


(73) 


which  when  combined  with  the  fundamental  dimensions  from  Table  3,  gives 

n  =  (ml3  f  (j L)a 2  (Z,)“3  (lt~{  f  '  (it1  Y  (milt  2  Y  (ml2t~2  Y 
(ml  [t~2  Y  (ml1t  2  Y  (ml~1t  2  Y°  {t )a"  '  (74) 

This  can  be  simplified  to 

n  =  (m)a  {l)Pi  (r)A 

where, 

/3\  —  Of]  +  CL  6  +  CC-j  +  CX 8  +  OCq  +  OfjQ ,  (75 ) 

P2  =  -3 al  +  a 2  +  a2  +  a4  +  a5  +  2 a6  +  2 an  -  ag  -  a9  -  al0 , 

J33  =  -a4  -  a5  -  2 ab  -  2 an  -  2 a&  -  2 a9  -  2al0  +  an. 


Setting  the  values  of  [3 to  zero  and  solving  for  ai,  0C2,  and  ogives 


ax  =  —  a6  —  cc-f  —  ctg  —  a9  —  al0 , 

a4  =  - a5  - 2 a6  - 2 a7  - 2 a8  - 2 a9  - 2a l0  +an, 

a,  =  —  a3  —  3  a6  —  3  an  —an. 


(76) 


This  set  of  equations  allows  the  definition  of  the  invariants.  Rewriting  these  equations  in 
vector  form  gives 
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Each  of  these  columns  represents  a  separate  invariant  that  must  be  met.  The  invariant  is 
found  by  associating  each  dimensioned  quantity  with  its  corresponding  a  value  and 
raising  the  dimensioned  quantity  to  the  power  seen  in  the  column  vector.  In  this  case,  the 
invariants  are  given  by 


Since  this  must  be  invariant  regardless  of  the  arbitrary  values  of  C,  the  separate  invariants 
are  found  by  setting  one  C  to  one  and  the  others  to  zero,  which  gives  the  invariants  as. 
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Kx 


d 


K1 


C 

V 

E 

2/3 

pv  d 
S 

pv2d 3 


<7 


^5  =' 


7.c 


^6  =' 


7T~]  -  ' 


pv 

K 

pv2 

Go_ 

pv2 


tv 

71  a  - 

8  d 


(79) 


This  process  has  now  given  a  set  of  physical  invariants  that  must  be  satisfied  in  every 
case  to  ensure  that  a  scaled  model  of  a  cylindrical  impact  is  equivalent  to  the  initial 

model.  The  last  invariant,  tt8  =  —  ,  is  used  in  comparing  two  different  scaled  models. 

d 

By  solving  this  invariant  for  the  time,  it  is  possible  to  detennine  which  time  step  in  model 
one  is  equivalent  to  which  time  step  in  model  2. 


3.3  Model  Scaling 

A  typical  test  sled  at  the  Holloman  High  Speed  Test  Track  consists  of  a  rocket 
attached  to  a  test  sled,  which  in  turn  is  attached  to  the  narrow  gauge  rail  via  either  four  or 
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six  shoes,  as  was  seen  in  Figure  1  and  Figure  2.  An  entire  sled  system  may  weigh  as 
much  as  809  kg.  The  shoes  that  connect  the  sled  to  the  rails  are  generally  20.32  cm  long, 
by  10.8  cm  wide,  by  2.54  cm  high.  The  Szmerekovsky  model  uses  a  full  size  shoe,  and 
models  the  mass  of  the  sled  as  being  evenly  distributed  across  four  shoes.  Because  CTH 
is  an  Eulerian  model,  and  not  a  Lagrangian  model,  it  is  not  possible  to  add  the  mass  as  a 
point  mass  at  a  node.  This  led  to  the  addition  of  an  artificial  sled  mass  on  top  of  the  shoe 
to  represent  the  quarter  mass  of  the  sled.  The  last  step  in  developing  the  Szmerekovsky 
model  involved  taking  a  unit  width  slice,  to  represent  a  plane  strain  model.  The  flow 
chart  showing  how  he  simplified  the  system  is  seen  in  Figure  24. [2] 


Figure  24  -  Sled  Mass  Simplification 

To  mitigate  gouging,  it  is  desirable  to  know  exactly  what  the  conditions  are  at  the 
moment  a  gouge  forms.  Ideally,  the  best  way  to  perform  this  would  be  to  record 
conditions  as  a  gouge  occurs,  on  the  actual  test  track.  The  problem  with  doing  this  is  that 
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high  velocity  impact  experiments  on  the  HHSTT  are  very  expensive  to  perform  in  terms 
of  cost,  time,  and  logistics.  In  addition  to  these  difficulties,  it  is  not  currently  possible  to 
determine  exactly  where  a  gouge  will  develop,  and  therefore  it  would  not  be  necessary  to 
instrument  the  whole  track  and  record  the  conditions  that  occur  during  a  gouging  event, 
which  is  not  a  viable  option.  Due  to  these  facts,  it  is  not  feasible  to  study  the  gouging 
phenomena  using  the  actual  test  track.  Therefore,  it  is  desirable  to  study  the  phenomena 
either  numerically,  or  using  an  equivalent  experiment  which  approximates  the  conditions 
seen  in  the  field.  Laird  and  Szmerekovsky  have  developed  a  numerical  model  that  scales 
the  physical  problem  down  to  a  level  that  it  can  be  simulated  using  a  reasonable  cluster  of 
Linux  computers.  [1,2] 

The  Szmerekovsky  model  did  an  excellent  job  of  simulating  a  gouging  event,  and 
provided  excellent  insight  into  the  formation  of  a  gouge  on  a  realistic  model.  Gouges 
modeled  using  this  model  were  found  to  be  similar  in  shape  and  size  to  those  seen  at  the 
HHSTT.  The  problem  faced  by  Szmerekovsky  however  was  that  the  physical  gouges 
seen  at  the  HHSTT  were  not  accompanied  by  any  description  of  the  conditions  under 
which  the  gouge  formed.  Without  information  on  the  conditions  of  the  impact,  it  was 
impossible  to  know  whether  the  impact  modeled  was  similar  to  the  impact  that  caused  the 
physical  gouge  seen  in  the  rail  at  Holloman. 

This  research  attempts  to  develop  a  simplified  gouging  model  that  is  equivalent  to 
the  Szmerekovsky  model,  and  therefore  the  physical  test  sled,  yet  is  simple  enough  to 
verify  using  an  oblique  ballistic  impact  experiment.  Most  high  velocity  impact 
experiments  are  performed  using  a  compressed  gas  gun  that  shoots  either  a  long 
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cylindrical  rod  projectile,  or  a  spherical  projectile.  As  will  be  seen  shortly,  it  is  important 
to  match  diameter  to  length  ratios  of  models.  Since  the  Szmerekovsky  model  utilizes  a 
rectangular  slipper  with  an  aspect  ration  of  1/8,  a  spherical  projectile  is  not  possible,  and 
a  rectangular  impact  rod  must  be  used.  It  was  further  decided  that  simply  using  a  right 
circular  cylinder  would  introduce  too  great  of  a  discontinuity  when  the  90  degree  corner 
of  the  cylinder  struck  the  impact  plate.  Because  of  this,  it  was  decided  to  use  a  spherical 
tip  on  the  impacting  end  of  the  rod. 

Szmerekovsky ’s  work  presented  a  number  of  different  impact  models  used  to 
represent  the  formation  of  a  gouge.  Models  were  developed  for  an  oblique  impact,  for  a 
horizontal  impact  with  a  circular  asperity  on  the  rail,  and  for  the  case  of  an  impact  with 
an  elliptical  “rail  roughness”,  as  may  be  seen  on  the  test  track.  Additional  cases  were 
also  investigated  using  an  epoxy  coating,  which  is  similar  to  the  coating  used  on  the  rails 
at  the  HHSTT.  Since  this  investigation  was  to  develop  a  model  that  could  be 
experimentally  verified,  it  was  determined  that  both  an  asperity  impact  and  a  rail 
roughness  impact  were  too  difficult  to  perform  in  a  controlled  manner  in  a  laboratory. 
Additionally,  it  was  decided  to  investigate  only  the  case  of  an  uncoated  rail. 

A  Buckingham  Pi  analysis  was  performed  based  upon  Szmerekovsky ’s  sled 
system  mass  model.  This  model  consisted  of  a  full  sized  slipper  with  an  artificial  mass 
added  on  top  to  simulate  the  mass  of  the  sled,  and  was  compared  with  a  long  rod  impact 
specimen  capable  of  being  launched  by  a  typical  compressed  gas  gun.  Szmerekovsky’ s 
model  and  the  model  developed  here  are  seen  in  Figure  25. [2] 
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Figure  25  -  Impact  Model  Comparison 

The  analysis  was  based  upon  Szmerekovsky’s  conservation  equation  case.  The 
dimensioned  quantities  used  are  given  in  Table  4.  From  this  it  is  seen  that  the 
fundamental  dimensions  are  mass  M,  length  L,  and  time  T. 

The  invariant  parameter  II  is  found  to  be, 

n  =  (mf  iff  (df  (wf  («,)"’  (uy  Y  (cf  (Bp  (Sp  (<v  f“  (E.  P  (G„  f  (tf  (SO) 
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Table  4  -  Buckingham  Pi  Dimensioned  Quantities  for  Rod  Impact  Model 


Dimensioned  Quantity 

Symbol 

Fundamental 

Dimensions 

Szmerekovsky 

Rod  Impact 
Model 

Szmerekovsky 

Rod 

Impact 

Model 

height  of 
slipper 

diameter  of  rod 

h 

d 

L 

length  of 
slipper 

length  of  rod 

/ 

1 

L 

horizontal 

velocity 

horizontal 

velocity 

ux 

Ux 

LT'1 

vertical 

velocity 

vertical 

velocity 

Uy 

Uy 

LT'1 

material  speed 
of  sound 

material  speed 
of  sound 

C 

c 

LT"1 

internal  energy 

internal  energy 

E 

E 

ML2T'2 

energy  source 

energy  source 

S 

S 

ML2T'2 

yield  strength 
in  compression 

yield  strength 
in  compression 

Gy,c 

Gy,c 

ml't2 

elastic 

modulus 

elastic  modulus 

Em 

f 

ML'T2 

shear  modulus 

shear  modulus 

G0 

G0 

ml't2 

time 

time 

t 

t 

T 

This  can  be  simplified  to 


n  =  (m)a  ( lY 2  ( t)p 3 

where, 

fi\  ~  &\  T  ^8  ^ 9  aX0  ^11  ^12  ’  (82 ) 

fi2  =  a2  +  a3  +  a4+  a5  +a6+a2  +  2 ag  +  2 a9  -  al0  -  an  -  aX2 , 
fi2  =  -a5  -a6— a -j  —  2  a%  -  2  a9  -2aw  -2an  -2al2  +al3. 


Setting  the  values  of  (3  to  zero  and  solving  for  ay,  a,2,  and  a-,- gives 
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(83) 


=  —  a8  —cc9  —a  10  —  ccn  —ccn, 
a5  =  —a6  -a7  -2 ai  -2 a9  -2 a10  —2 an  -2 au  +a13, 
a2  =  -a3  - a4  +3 aw  +3 au  +3a12  -an. 

This  set  of  equations  allows  the  definition  of  the  invariants.  Rewriting  these  equations  in 
vector  form  gives 
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Each  of  these  columns  represents  a  separate  invariant  that  must  be  met.  The  invariant  is 
found  by  associating  each  dimensioned  quantity  with  its  corresponding  a  value  and 
raising  the  dimensioned  quantity  to  the  power  seen  in  the  column  vector.  In  this  case,  the 
invariants  are  given  by 


n  = 


w  j 


f-T' 

l/J 


(  _  A 


vu,y 


Kmux  J 


Q  / 


^9/ 


V  mUx  J 


a  I 

y,c 


3  ^'°'E  l^C" 


V  mUx  J 


V  mUx  J 


ymux  j 


(85) 


\  1  J 


Since  this  must  be  invariant  regardless  of  the  arbitrary  values  of  C,  the  separate  invariants 
are  found  by  setting  one  C  to  one  and  the  others  to  zero,  which  gives  the  invariants  as. 
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nx 

n2 


7T  ^ 


ka 


X5 


7V 


6 


n7 


7X% 


7T  g 


n 


10 


d 

/ 

w 

T 

"y_ 

UX 

c 

u 


mu . 

S 

2 

mux 

°yf 

2 

mUx 

E  /3 

m 

2 

mUx 

Ga 

2 

mUx 

tu^ 

1 


(86) 


To  maintain  proper  scaling  these  invariant  Pi  parameters  must  be  matched  in  both 
the  Szmerekovsky  model  and  the  rod  impact  model.  The  values  of  the  dimensioned 
quantities  as  originally  used  in  the  Szmerekovsky  model  are  shown  in  Table  5. 
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Table  5  -  Szmerekovsky’s  Dimensioned  Quantities 


mass 

19.1  kg 

length 

20.32  cm 

height 

2.54  cm 

width 

10.8  cm 

horizontal  velocity 

3000  m/s 

vertical  velocity 

-2  m/s 

material  speed  of  sound 

4737  m/s 

internal  energy 

930  Mega  Joules 

energy  source 

0  Mega  Joules 

yield  stress  in  compression 

14.47  GPa 

elastic  modulus 

1824  GPa 

shear  modulus 

718  GPa 

Analysis  was  performed  using  a  two-dimensional  plane  strain  solution,  which  called  for 
the  use  of  a  unit  width  of  the  slipper.  Since  a  unit  width  was  used,  the  total  mass  of  the 
system  was  divided  by  the  width  of  the  slipper.  The  total  mass  is  found  by  adding  one 
quarter  of  the  mass  of  the  sled  to  the  mass  of  the  slipper,  which  is  found  by  calculating 
the  volume  of  the  slipper  and  multiplying  by  the  density  of  VascoMax  300,  which  is 
8.129  g/cm3.  The  material  speed  of  sound  is  calculated  from  the  equation  presented  in 
[E~ 

Chapter  2,  c  =  .  Prior  to  the  impact,  the  internal  energy  of  the  system  is  entirely  in 

V  P 

1 

the  kinetic  energy  of  the  projectile,  which  is  —  mv  ■  There  is  no  internal  energy  source 

within  the  projectile,  so  therefore  S  is  zero. 

The  values  calculated  for  the  Pi  parameters  from  the  Szmerekovsky  model  are, 
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kx  =  .125 
n2  =  .531 
tt3  = -6.667  xlO4 
nA  =3.158 

tt5  =0.500  .  (87) 

x6=0 

7T1  =  .282 

n%  =35.519 

;r9  =  13.984 

The  model  to  be  investigated  must  now  be  developed.  In  order  to  most  closely 
meet  the  goals  of  studying  the  gouging  phenomena,  it  was  desirable  to  use  the  actual 
materials  whenever  possible.  Because  of  this,  it  was  assumed  that  the  projectile  was 
made  of  VascoMax  300  steel  with  the  same  properties  as  in  the  Szmerekovsky  model. 
This  constrained  the  yield  strength,  elastic  modulus,  and  shear  modulus,  as  well  as  the 
density,  even  though  this  is  not  used  directly  in  the  Buckingham  Pi  process.  With  the 
material  properties  set,  it  was  left  to  determine  the  dimensions  of  the  projectile  and  the 
velocity  of  impact.  Since  the  original  materials  are  used,  the  sound  speed  is  constant 
between  the  two  models,  while  the  horizontal  velocity  is  varied.  Because  of  this,  it  is 
necessary  to  allow  714  to  vary  between  the  two  models. 

Most  impact  projectiles  used  are  fairly  small,  because  the  larger  the  projectile,  the 
greater  the  amount  of  energy  that  is  needed  to  accelerate  it  to  speed.  In  general,  most 
projectiles  are  between  5  and  10  mm  in  diameter.  Because  of  this,  the  diameter  of  the  rod 
was  initially  chosen  to  be  6  mm.  Starting  from  this  diameter,  and  using  Tli  from 
Szmerekovsky ’s  model  gives  a  rod  length  of  48  mm.  The  volume  of  the  model  of  the 
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rod,  which  is  also  a  plane  strain  model  with  a  unit  width,  was  found  by  adding  a 
rectangular  rod  and  a  half  circle  for  a  tip,  and  multiplying  by  a  unit  width, 


volume  = 


,  ,  1 

fd) 

2 

a  •/  +  —  n 

2 

U ) 

(88) 


Using  the  volume  and  the  density,  the  mass  was  calculated,  which  left  the  velocity  to  be 
determined. 

The  velocity  was  determined  using  a  scaling  factor.  To  begin,  an  initial  guess  of 
1000  m/s  was  chosen  for  the  horizontal  velocity  ux.  From  this  initial  guess  of  velocity, 
the  invariant  n8  was  calculated  using 

/3 

8  )  GUESS'  —Em  /  y  2  (8$ 

(p  ■  volume  )ux 


where  the  term  in  parentheses  is  the  mass  of  the  rod.  This  result  was  then  compared  to 
the  original  II8  as  calculated  from  Szmerekovsky’s  model,  and  a  ratio  was  developed  that 
was  the  scaling  factor, 


scale  = 


GUESS 

n8 


(90) 


This  scaling  factor  is  utilized  along  with  the  definition  of  n8  to  detennine  the  actual 
horizontal  velocity  of  the  projectile. 
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scale  = 


(91) 


_ /3 _ 

p*  volume  *(uxGUESS)2 

_ /3 _ 

p*  volume *(ux  ACTUAL)2 


Because  both  Pi  parameters  are  based  upon  the  same  model,  with  different  velocities,  it  is 
possible  to  simplify  this  equation  and  solve  for  uXiactual  which  gives  the  horizontal 
velocity  that  must  be  used, 

ux, actual  ~  v/ scale  •  ux  GUESS .  (92) 


Determination  of  the  horizontal  velocity  allows  the  calculation  of  the  rest  of  the  Pi 
parameters,  which  are  given  in  Table  6,  and  compared  to  those  of  Szmerekovsky’s 
model.  Additionally,  Pi  parameters  are  given  for  both  the  Szmerekovsky  model  and  the 
rod  impact  model  for  Szmerekovsky’s  case  of  a  horizontal  velocity  of  1500  m/s  and 
vertical  velocity  of  1  m/s. 


Table  6  -  Pi  Parameters  for  Szmerekovsky  Model  and  Rod  Model 


ux=3000  m/s 
uv=-2  m/s 

ux=1500  m/s 
uv=-l  m/s 

Szmerekovsky 

Rod  Model 

Szmerekovsky 

Rod  Model 

Ux 

3000  m/s 

9617  m/s 

1500  m/s 

4808.5  m/s 

Uv 

-2  m/s 

-6.41  m/s 

-1  m/s 

-3.21  m/s 

Th 

0.125 

0.125 

0.125 

0.125 

n3 

-6.667x1 0'4 

-6.667xl0'4 

-6.667x1 0'4 

-6.667x1 0'4 

n4 

1.579 

0.493 

3.158 

0.985 

n5 

0.5000 

0.5000 

0.5000 

0.5000 

n6 

0 

0 

0 

0 

n7 

0.07045 

0.07045 

0.282 

0.282 

n8 

8.88 

8.88 

35.519 

35.519 

n9 

3.496 

3.496 

13.984 

13.984 
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In  order  to  keep  materials  constant  between  the  two  models,  it  was  necessary  to  allow  the 
invariant  relating  material  sound  speed  and  horizontal  velocity  to  be  non-constant.  Any 
attempt  to  perfectly  constrain  all  invariants  would  have  resulted  in  a  change  of  material 
properties  that  would  have  been  impossible  to  duplicate  in  the  field. 


Chapter  4  -  Analysis  and  Results 


4.1  Verification  of  CTH  Model  Parameters 

Prior  to  embarking  on  an  effort  to  model  an  oblique  impact,  it  was  necessary  to 
determine  the  proper  parameters  to  use  in  CTH.  To  do  this,  a  model  was  developed 
based  upon  the  Taylor  cylinder  test.  The  Taylor  cylinder  test  involves  striking  a  rigid 
anvil  at  a  right  angle  with  a  right  circular  cylinder,  as  was  seen  in  Figure  7.  This  test  is 
generally  used  to  determine  the  dynamic  yield  stress,  [10]  however  in  this  case,  the 
pictures  obtained  from  the  work  of  Jones,  et  ah,  [14]  are  compared  to  the  final  deformed 
specimen  as  simulated  using  a  Taylor  cylinder  test  model  built  for  CTH. 

The  model  built  for  CTH  was  based  upon  the  experiment  performed  by  Jones,  et 
al.  In  their  experiment,  Jones  used  an  Oxygen  Free  High  Conductivity  (OFHC)  copper 
impact  cylinder  with  a  diameter  of  7.62  mm  and  a  length  of  57. 15  mm.  This  cylinder  was 
fired  at  a  velocity  of  176  m/s  against  a  4340  steel  anvil.  Fortunately,  the  materials  library 
in  CTH  contains  the  equation  of  state  and  the  constitutive  equation  for  both  copper  and 


4340  steel;  therefore  these  materials  could  be  modeled  perfectly  in  the  simulation.  [14] 
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The  boundary  condition  for  the  edge  of  the  mesh  that  was  used  was  the 
hydrodynamic  boundary  condition.  Under  this  condition,  when  a  stress  wave  reaches  a 
mesh  boundary,  it  is  imagined  to  continue  unimpeded  as  if  the  material  extended 
infinitely.  The  only  instance  where  stress  waves  reflect  is  when  there  is  a  change  in 
materials,  such  as  at  the  material  interface  between  the  rod  and  plate,  or  along  the  sides  of 
the  rod  where  the  rod  material  ends  and  a  void  is  encountered.  The  mesh  used  in  the 
Taylor  test  simulation  was  0.025  cm  square  in  the  region  of  impact,  and  throughout  the 
rod.  In  the  area  away  from  the  impact,  a  mesh  of  0. 100  cm  was  used.  A  diagram  of  the 
mesh  is  presented  in  Figure  26.  Only  half  of  the  Taylor  test  needs  to  be  modeled, 
because  it  can  be  represented  axisymetrically,  and  the  other  half  of  the  model  can  be 
added  using  a  mirror  option.  The  target  was  chosen  as  a  20  cm  thick  piece  of  1080  steel. 
This  is  much  more  than  is  needed,  since  the  boundary  conditions  were  essentially  semi- 
infinite,  and  therefore  the  stress  wave  would  not  be  reflected  no  matter  how  long  the 
simulation  was  carried  out. 
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x=0.381  cm 


Figure  26  -  Taylor  Test  Grid 

It  was  left  to  be  determined  which  user  controllable  options  within  CTH  should  be 
utilized.  The  major  options  within  CTH  which  can  be  changed,  and  may  change  results, 
are  the  method  for  determining  yield  strength  in  mixed  cells,  the  handling  of  multiple 
materials  and  pressures,  whether  voids  within  the  specimen  should  be  given  any  strength, 
and  which  material  interface  algorithm  should  be  used,  the  algorithm  that  moves  the 
sliding  away  from  the  material  interface,  or  the  algorithm  that  sets  the  deviatoric  stress  to 
zero. 

A  baseline  model  was  developed  using  the  following  options,  which  will  be 
briefly  described: 

•  Yield  strength  in  mixed  cells  -  The  yield  strength  in  a  cell  containing 
multiple  materials  is  the  volume  averaged  yield  strength  of  all  materials  in 
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the  cell  normalized  by  the  sum  of  the  volume  fractions  of  materials  that 
support  shear. 

•  Material  volume  fractions  and  pressures  in  mixed  cells  -Multiple 
materials  in  a  given  cell  can  each  have  their  own  temperature  and 
pressure,  however  there  is  no  means  of  relaxing  the  pressure  within  the 
cell. 

•  Void  strength  -  With  this  option,  when  a  void  is  developed  in  a  material 
due  to  local  fracture,  the  void  is  compressed  closed  as  if  it  has  no  strength. 

•  Interface  layer  -  This  option  creates  a  pair  of  boundary  layers  at  the 
material  interface,  one  hard  and  one  soft.  In  essence,  the  sliding  that 
would  occur  at  the  interface  of  the  two  materials  is  moved  into  the  soft 
boundary  layer,  effectively  moving  the  sliding  interface  of  the  two 
materials  into  the  “soft”  material. 

There  are  additional  options  that  could  have  been  investigated,  however  it  was  thought 
that  those  additional  options  were  well  enough  understood  to  choose  correctly  without 
investigation. 

The  CTH  code  was  used  to  simulate  the  Taylor  Test  perfonned  by  Jones. [14]  The 
final  deformed  shape  calculated  by  CTH  was  compared  to  the  deformed  shape  as 
photographed.  Figure  27  shows  that  the  baseline  CTH  simulation  at  least  approximately 
simulates  the  experimental  results  near  the  impacted  end.  Using  the  baseline  model  as  a 
starting  point,  each  of  the  four  input  parameters  that  could  be  controlled  were  changed 
one  at  a  time  to  detennine  which  options  to  use. 
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(a)  (b) 

Figure  27  -  Taylor  Test  Specimen  -  (a)  Experimental  (b)  CTH  baseline 


The  first  option  that  was  adjusted  was  the  calculation  of  the  yield  strength  in 
mixed  cells.  Initially,  the  yield  strength  was  calculated  based  upon  the  volume  averaged 
yield  strength,  normalized  by  the  sum  of  the  volume  fractions  that  can  support  shear.  The 
other  option  investigated  sets  the  strength  in  mixed  cells  to  zero,  except  if  a  cell  contains 
one  material  and  a  void.  It  would  be  suspected  that  this  might  give  unrealistically  weak 
material  response  near  the  interface  where  mixing  was  expected.  The  CTH  simulation 
showed  that  this  appeared  to  be  the  case,  as  even  at  the  low  speeds  of  this  simulation, 
particles  of  material  were  ejected  from  the  specimen,  as  can  be  seen  in  Figure  28.  This 
option  did  not  improve  the  model  and  is  not  used  further. 
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Material  Impact 


-2 . 0 

•-••• ,  -2  -1  0  1  2 

2DC  X  (cm) 

Taylor  Test:  Capper  on  4340  Steel,  V=176  m/s,  6cm  specimen 

Figure  28  -  Taylor  CTH  Simulation  -  Yield  Strength 
The  manner  in  which  CTH  calculates  pressure  in  cells  of  mixed  materials  was 

also  investigated.  The  baseline  model  allowed  materials  to  have  independent  pressures 

within  a  mixed  cell,  however  it  doesn’t  allow  pressure  relaxation.  The  most  sophisticated 

of  the  multiple  material  and  pressure  models  allows  multiple  materials  and  pressures  in 

mixed  cells,  proportional  to  the  material  volume  fraction,  however  it  does  not  allow 

volume  change  for  materials  with  small  volume  fractions.  This  allows  pressure  to  relax 

to  equilibrium  values.  Since  this  is  a  more  sophisticated  model,  it  was  expected  that 

results  would  be  as  good  or  better  than  the  baseline  model.  Figure  29  presents  the 

deformed  shape  of  the  Taylor  specimen  using  the  new  option.  Qualitative  observation 

shows  that  this  model  is  superior  to  the  previous  model,  because  it  seems  to  allow  the 

material  farthest  from  the  rod’s  axis  to  flow  more  freely  as  it  appears  to  do  in  the 


73 


experimental  test.  This  flow  of  material  results  in  a  sharper  edge  on  the  mushroomed 
head  than  appears  in  the  baseline  model,  which  appears  to  be  slightly  blunted. 


Material  Impact 


Taylor  Test:  Copper  on  4340  Steel,  V=176  m/s,  6cm  specimen 

Figure  29  -  Taylor  CTH  Simulation  -Multiple  Materials  and  Pressures 
The  baseline  CTH  model  allowed  voids  that  form  within  the  specimen  to  be 

compressed  from  the  specimen  before  the  surrounding  material  was  compressed.  An 

option  exists  which  in  effect  gives  the  void  within  the  material  a  strength,  because  it 

compresses  voids  and  material  according  to  their  volume  fractions.  This  requires  energy 

to  deform  the  material,  which  might  otherwise  be  used  in  closing  voids.  Results  for  this 

option  also  ejected  material  from  the  projectile,  as  seen  in  Figure  30.  Due  to  this  fact, 

voids  will  be  allowed  to  close  as  in  the  baseline  model. 
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Taylor  Test:  Copper  on  4540  Steel,  V=176  m/s,  6cm  specimen 


Figure  30  -  Taylor  CTH  Simulation  -  Void  Strength 
The  last  option  investigated  was  the  manner  in  which  the  boundary  between  the 


two  materials  was  simulated.  The  baseline  CTH  model  utilized  the  algorithm  that 


effectively  moves  the  sliding  interface  into  the  softer  material.  The  other  option 


investigated  sets  the  shear  force  at  the  interface  to  zero.  The  second  algorithm  appears  to 
give  good  results  as  well,  as  can  be  seen  in  Figure  3 1 . 
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Material  ImDact 


Taylor  Test:  Copper  on  4340  Steel,  V=176  m/s,  6cm  specimen 


Figure  31  -  Taylor  CTH  Simulation  -  Interface  Layer 
The  appropriate  options  to  use  for  a  Taylor  Test  simulation  were  just  shown  to  be 

allowing  multiple  materials  and  pressures  in  a  cell,  with  pressure  relaxation,  and  also  to 

set  the  shear  stress  at  the  boundary  to  zero.  These  two  changes  were  made  to  the  baseline 

model,  and  an  additional  CTH  simulation  was  perfonned.  This  model  was  compared  to 

the  Jones  experiment  through  the  use  of  tracer  points  in  the  CTH  simulation.  Tracer 

points  are  Lagrangian  points  within  the  mesh  that  are  followed  throughout  deformation, 

allowing  the  material  history  of  a  specific  point  to  be  followed.  For  the  simulation,  sixty 

tracer  points  were  placed  on  the  free  edge  of  the  cylinder,  along  the  axis.  After 

deformation,  the  positions  of  these  points  were  graphed,  presenting  the  exact  deformation 

of  the  rod.  This  deformed  shape  was  compared  to  the  measured  profile,  as  presented  in 

the  paper  by  Jones  et  al.[  10]  The  deformed  shapes  of  the  experimental  and  numerical 

methods  are  presented  in  Figure  32.  For  simplicity,  only  one  half  of  the  rod  is  shown, 
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and  the  scales  for  radius  and  length  are  not  the  same,  which  is  why  the  rod  appears 
squatter  than  in  the  picture  presented  earlier.  The  defonned  shape  seen  in  Figure  32, 
shows  that  the  use  of  multiple  materials  and  pressure  algorithm  and  interface  algorithm 
previously  described,  is  justified,  because  the  simulated  shape  is  very  similar  to  the  actual 
deformed  shape  recovered  from  an  experimental  setup.  The  axial  position  where  both 
tests  begin  to  mushroom  appears  to  be  approximately  25  mm  from  the  impacted  end.  The 
lip  seen  in  the  experimental  test  at  around  3  mm  from  the  impacted  end  is  simulated 
nearly  exactly  in  the  CTH  simulation.  Pictures  of  the  experimental  impact  specimen  and 
the  CTH  simulated  specimen  are  seen  in  Figure  33. 

The  two  major  areas  where  the  two  experiments  differ  are  at  the  two  ends.  The 
CTH  model  over  estimates  the  total  amount  of  axial  defonnation  by  approximately  1.3 
mm,  which  is  only  3%  of  the  length  of  the  deformed  specimen.  The  other  are  where  the 
two  tests  differ  is  at  the  impacted  end,  where  the  CTH  model  underestimates  the 
deformation  by  approximately  1  mm.  This  difference  can  be  explained  by  the  fact  that 
the  interface  algorithm  does  not  permit  any  friction  at  the  interface.  This  would  seem  to 
cause  the  CTH  model  to  overestimate  the  deformation,  however  since  there  is  no  friction 
in  this  region,  any  temperature  change  is  due  solely  to  plastic  deformation.  Since  there  is 
no  deviatoric  stress  immediately  at  the  interface,  the  temperature  rise  is  due  solely  to 
plastic  defonnation,  and  will  likely  be  much  lower  than  in  experiments.  The  lower 
temperature  will  result  in  less  plastic  flow,  because  the  yield  stress  won’t  decrease  as 
much  in  the  constitutive  equation,  and  therefore  the  yield  stress  won’t  be  reached  as 
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early,  and  less  plastic  flow  will  occur.  With  less  plastic  flow  occurring,  the  diameter  of 
the  impacted  end  will  tend  to  be  smaller  than  in  the  experimental  case. 


Deformed  Rod 


Figure  32  -  Taylor  Test  -  CTH  Model  Versus  Experimental  Results 
Qualitative  evaluation  of  the  CTH  model  presented  over  the  previous  few  pages 

shows  that  CTH  generates  a  deformed  Taylor  Test  specimen  that  is  very  similar  to  an 

experimental  Taylor  Test  specimen.  It  was  left  to  determine  if  CTH  was  quantitatively 

accurate  as  well.  The  theory  discussed  in  Section  2.2  describes  some  of  the  behaviors 

that  should  be  witnessed  using  the  CTH  simulation. 
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(a)  (b) 

Figure  33  -  Deformed  Taylor  Test  (a)  Experiment  (b)  CTH 


The  one-dimensional  theory  of  stress  in  a  long  rod  subjected  to  an  impact,  is  an 
approximate  theory  that  states  that  after  impact,  a  one-dimensional  stress  wave  is  formed 

jp 

that  travels  at  constant  speed  equal  to  c  =  I-2-  where  Em  is  the  elastic  modulus  of  the 

V  P 

material  with  density  p.  The  stress  state  behind  the  wave  is  one  of  uniaxial  compression 
where  the  stress  is  constant  and  equal  to  cr  =  pcvo  where  vQ  is  the  impact  velocity.  These 
theories  were  checked  utilizing  a  Taylor  Test  impact  of  high  strength  VascoMax  300 
steel  on  1080  steel,  at  a  velocity  of  176  m/s.  The  stress  profde  at  0.5  ps  is  seen  in  Figure 
34.  This  picture  shows  the  formation  of  the  one-dimensional  stress  wave  predicted  by 
theory.  Elementary  one-dimensional  impact  theory  predicts  that  the  compressive  stress 
will  be  770  MPa  behind  the  stress  wave.  Figure  34  shows  that  the  stress  wave  is 
approximately  700  MPa. 
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Taylor  Test:  V300  on  1080  steel,  v=176  m/s,  6cm  specimen 
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Figure  34  -  Stress  Wave  at  0.5  jus 

The  stress  wave  is  assumed  to  travel  at  a  constant  speed  of  4,930  m/s  through  the 


steel  rod,  which  is  6.0  cm  long,  which  means  that  it  should  take  12.4  ps  to  travel  the 


length  of  the  rod.  The  rod  is  shown  at  12.5  ps  in  Figure  35.  If  the  stress  wave  is  assumed 
to  reach  the  end  of  the  rod  shortly  after  it  is  seen  in  this  figure,  it  can  be  assumed  that  it 
took  approximately  12.5  ps  to  traverse  the  rod,  which  is  as  expected.  It  is  difficult  to 
exactly  determine  when  the  stress  wave  reaches  the  end  of  the  rod  due  to  the  time  step 
utilized,  and  the  fairly  coarse  discretization  of  the  rod  far  from  the  impacted  end.  The 
reason  that  the  stress  farther  from  the  wave  front  is  greater  than  700  MPa  is  that  in  the 
theory  of  long  rod  impact,  it  is  assumed  that  the  stress  doesn’t  interact  with  the  sides  of 
the  rod,  when  in  fact  the  stress  will  reflect  from  the  sides,  changing  the  stress  state  from 
one  of  strictly  one-dimension  to  a  state  of  three-dimensional  stress. 
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Taylor  Test:  V300  on  1080  steel,  v=176  m/s,  6cm  specimen 


Figure  35  -  Stress  Wave  at  12.5  /us 

The  last  piece  of  theory  discussed  in  regards  to  stress  wave  propagation  in  long 
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rods  is  that  a  stress  wave  reflects  with  the  opposite  sign,  effectively  zeroing  the  stress 


behind  the  reflected  stress  wave.  This  zeroing  of  stress  can  clearly  be  seen  in  Figure  36, 
where  the  white  stress  contours  at  the  top  of  the  rod  are  approximately  zero  stress,  behind 


the  stress  wave,  which  is  moving  down. 
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Taylor  Test:  V300  or  1080  steel,  v=176  m/s,  6cm  specimen 


Figure  36  -  Stress  Wave  at  17.5  jus 


4.2  Oblique  Rod  Impact  Model  1500  m/s  and  3000  m/s 

The  oblique  rod  impact  model  described  in  Section  3.3  was  utilized  in  simulations 
of  a  3000  m/s  horizontal  by  2  m/s  vertical  sled  run  and  a  1500  m/s  horizontal  by  1  m/s 
vertical  sled  run.  The  two  velocities  were  used  to  match  the  test  cases  that  Szmerekovsky 
used,  which  were  chosen  because  in  general  the  HHSTT  begins  to  see  gouging  around 
1500  m/s,  and  it  is  desired  to  eventually  perform  tests  at  the  HHSTT  at  speeds  of  3000 
m/s.  [2]  Simulations  were  performed  on  a  cluster  of  64-bit  Linux  computers,  using  either 
ten  or  twelve  processors.  On  average,  the  simulations  took  40  minutes  to  perform. 
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The  boundary  condition  used  in  the  rod  impact  model  was  the  hydrodynamic 
boundary  condition.  This  boundary  condition  treats  mesh  boundaries  as  being  semi¬ 
infinite,  and  therefore  a  stress  wave  that  reaches  a  boundary  continues  on  its  present 
vector,  without  reflecting.  Material  boundaries  are  treated  realistically,  and  waves  reflect 
as  dictated  in  theory.  In  the  area  of  the  impact,  the  mesh  was  0.0050  cm  square.  Farther 
from  the  impact  zone,  the  mesh  is  coarser.  The  number  of  cells  in  the  target  was  95,200, 
and  in  the  rod  there  were  28,1 10  cells.  The  mesh  is  shown  schematically  in  Figure  37, 
the  rod  is  shaded  dark  gray,  and  the  picture  is  not  to  scale. 


Figure  37  -  Rod  Impact  Model  Grid 


The  first  simulation  performed  was  the  simulation  of  the  1500  m/s  run.  The  rod 
impact  model  had  parameters  shown  earlier  in  Table  6.  The  rod  was  4.8  cm  long  by  0.6 
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cm  in  diameter,  fired  at  4808.5  m/s  horizontal  and  3.2  m/s  vertical.  For  this  case,  a  mesh 
slightly  larger  than  that  used  by  Szmerekovsky  was  utilized.  In  the  area  of  impact,  cells 
were  0.0050  cm  on  a  side,  as  opposed  to  the  0.0020  cm  that  Szmerekovsky  used. [2] 

Early  in  the  deformation  event,  voids  are  seen  in  the  1080  steel  layer,  beginning  around 
5.5  ps.  These  voids  appear  to  be  wear,  which  is  described  by  Bayer  as  “progressive 
damage  to  a  surface  caused  by  relative  motion  with  respect  to  another  substancc/'[l  5] 

An  example  of  this  apparent  wear  is  seen  in  Figure  38,  where  the  red  material  is  the 
VascoMax  300  steel,  the  yellow  material  is  the  1080  steel,  and  the  white  is  representative 
of  a  void. 
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Figure  38  -  Rod  Impact  Model  Wear  at  4. 8  km/s 
Szmerekovsky  and  Laird  [2,1]  described  a  hump  of  material  that  was  found  to 

build  as  a  gouge  begins  to  form.  This  is  the  precursor  that  causes  the  materials  to 

impinge  upon  each  other  and  begin  to  gouge.  The  rod  impact  model  developed  this 

characteristic  hump  of  material,  as  seen  in  Figure  39.  The  hump  continues  to  increase  in 

size  as  seen  in  Figure  40.  As  the  hump  continues  to  increase  in  size,  it  takes  the  classic 
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shape  seen  by  both  Szmerekovsky  and  Laird  in  their  work.  Szmerekovsky’s  gouge  is 
seen  in  Figure  41.  The  gouge  formed  using  the  rod  impact  model  is  seen  in  Figure  42. 

A  simulation  was  perfonned  with  a  finer  mesh  in  the  impact  region,  to  study 
whether  a  more  defined  mesh  would  improve  results.  The  mesh  size  in  the  region 
surrounding  the  impact  area  was  decreased  in  size  from  0.0050  cm  to  0.0025  cm.  The 
resulting  gouge  simulated  using  a  finer  mesh  is  seen  in  Figure  43,  but  does  not  show 


much  improvement  over  the  coarser  mesh. 
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Figure  39  -  Hump  formation  at  4.8  km/s 
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Figure  40-  Hump  Growing 
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Figure  41  -  Gouge  from  Szmerekovsky ’s  work 
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Figure  42  -  Rod  Impact  Model  Gouge  -  0. 0050  cm  mesh 
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Figure  43  -  Rod  Impact  Model  Gouge  -  0.0025  cm  mesh 
The  rod  impact  model  has  been  shown  to  develop  the  basic  shape  of  a  gouge  seen 
by  Szmerekovsky  and  Laird;  the  question  remained  as  to  whether  the  internal  properties 
that  were  present  during  the  gouging  event  were  similar  to  those  seen  by  Szmerekovsky. 
As  gouging  occurs,  the  deformation  in  the  vicinity  of  the  gouge  is  virtually  all  plastic 
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deformation.  The  plastic  strain  of  the  rail  found  by  Szmerekovsky  is  shown  in  Figure  44. 
It  is  seen  that  in  the  portion  of  the  rail  that  jets  into  the  slipper,  the  plastic  strain  is 
approximately  5.5.  This  very  closely  matches  the  plastic  strain  seen  using  the  rod  impact 
model,  where  the  plastic  strain  is  approximately  5.5,  as  seen  in  Figure  45.  Similarly,  the 
plastic  strain  in  the  rod  was  determined.  Szmerekovsky ’s  model  developed  a  plastic 
strain  of  approximately  5.5  in  the  slider,  as  can  be  seen  in  Figure  46.  The  rod  impact 
model  developed  a  plastic  strain  of  approximately  4.0,  as  is  shown  in  Figure  47.  The 
lower  plastic  strain  in  the  rod  compared  to  the  slider  used  by  Szmerekovsky,  may  be  due 
to  the  fact  that  there  is  much  less  energy  in  the  rod  impact  model  with  which  to  deform 
the  rod. 

The  time  scale  determined  in  the  Buckingham  Pi  analysis  of  the  Szmerekovsky 
model  and  the  rod  impact  model  was  13.5.  What  this  means  is  that  13.5  ps  in  the 
Szmerekovsky  model  is  equivalent  to  1.0  ps  in  the  rod  impact  model.  The  time  in  both 
Figure  44  and  Figure  45  is  approximately  6.0  ps.  At  this  point  in  the  Szmerekovsky 
model,  the  rod  impact  model  should  be  at  0.44  ps,  or  vice  versa,  at  this  point  in  the  rod 
impact  model,  the  Szmerekovsky  model  should  be  at  8 1 .0  ps.  This  shows  that  the  rod 
impact  model  does  not  exactly  match  the  time  scale  of  the  Szmerekovsky  model. 
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Figure  44  -  Plastic  Strain  (Rail)  -  Szmerekovsky  Model 
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Figure  45  -  Plastic  Strain  (Rail)  -  Rod  Impact  Model 
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Figure  46  -  Plastic  Strain  (Slider)  -  Szmerekovsky  Model 
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Figure  47  -  Plastic  Strain  (Rod)  -  Rod  Impact  Model 
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Szmerekovsky  showed  that  as  a  gouge  develops,  a  high-pressure  core  forms  at  the 
junction  of  the  two  materials  as  seen  in  Figure  48.  The  rod  impact  model  developed  a 
similar  high-pressure  core  in  the  same  general  area  as  in  the  Szmerekovsky  model,  see 
Figure  49.  Numerically  they  do  not  match  very  well,  which  is  understandable,  due  to  the 
small  set  of  variables  chosen  for  the  Buckingham  Pi  scaling  of  the  model.  With  more 
variables  involved,  solutions  would  improve  numerically.  Both  models  show  the  feature 
that  the  pressure  tends  to  be  centered  more  in  the  rail  material  (1080  steel). 

The  deviatoric  stress  component  tends  to  be  the  component  of  stress  that  is  the 
most  important  in  the  development  of  plasticity.  The  deviatoric  stress  component  was 
studied  by  Szmerekovsky,  and  the  stress  profde  was  as  seen  in  Figure  50.  The  deviatoric 
stress  component  in  the  rod  impact  model  was  shown  in  Figure  5 1 .  The  horizontal  stress 
contours  seen  in  the  Szmerekovsky  model  are  seen  in  the  rod  impact  model  as  well. 

These  horizontal  lines  are  likely  due  to  numerical  abnormalities  brought  about  by  the 
mesh  size.  The  rod  impact  model  shows  a  spike  in  the  deviatoric  stress  at  the  point  where 
gouging  is  initiated,  which  is  represented  by  the  orange  stress  contour. 

Similarly,  the  shear  stress  using  the  two  models  was  determined.  Szmerekovsky ’s 
model  is  seen  in  Figure  52,  and  the  rod  impact  model  is  shown  in  Figure  53.  The  two 
models  have  similar  profiles,  and  the  values  of  stress  are  quite  similar.  Both  models 
showed  an  area  of  high,  negative  shear  stress  in  the  rod  material  above  the  jet  of  the  rail 
material  into  the  rod.  This  high  negative  shear  stress  is  followed  by  a  high  positive  shear 
stress  farther  to  the  right  along  the  crest  of  the  jet. 
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Figure  48  -  High-Pressure  Core  -  Szmerekovsky 
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Figure  49  -  High-Pressure  Core  -  Rod  Impact  Model 


92 


XY  Deviators 


3  .  5 
3  .  0 
2  .  5 
2  .  0 


1  .  5 

>3-  1.0 

>- 

0  .  5 
0  .  0 
-0  .  5 

-1.0 


-1.5  L 
1  9 

2DR 


20  21  22 
X  (cm) 


Stress  Devlator  XY 
(d/ cm2) 

5.000X1C9 
4.524x10s 
4.048x10s 
3.571x10* 
3.095x10* 
2.619x10* 
2.143x10s 
1.867x10* 
1.190x10* 
7.143x10* 
2.381x10s 
1  -2.381x10“ 

I  -7.143x10s 
I  -1.190x10* 
-1.667x10* 
-2.143x10* 

I  -2.619x10* 


-3.095x10" 

-3.571x10* 

-4.048x10* 

-4.524x10* 

-5.000x10* 


2  3 


Figure  50  -  Deviatoric  Stress  -  Szmerekovsky 
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Figure  51  -  Deviatoric  Stress  -  Rod  Impact  Model  at  Gouge  Initiation 
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Figure  52  -  Shear  Stress  -  Szmerekovsky 
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Figure  53  -  Shear  Stress  -  Rod  Impact  Model 
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Szmerekovsky  also  investigated  the  motion  of  the  material  within  a  gouge.  He 
saw  that  in  the  vicinity  of  a  gouge,  the  material  tends  to  flow  in  the  general  direction  of 
the  gouge’s  formation.  His  gouge  model  is  shown  in  Figure  54,  with  the  velocity  vectors 
of  the  material  represented  by  arrows.  The  rod  impact  model  developed  a  very  similar 
velocity  profile,  as  seen  in  Figure  55.  The  rail  material  on  the  bottom  of  this  figure  has 
an  apparent  constant  velocity  to  the  left,  which  is  due  to  the  fact  that  the  mesh  is  given  a 
velocity  as  well,  to  maintain  gouging  in  the  finer  meshed  region 

Lastly,  the  temperature  was  investigated.  Szmerekovsky  found  a  temperature 
profile  as  seen  in  Figure  56.  The  temperature  profile  developed  using  the  rod  impact 
model  is  shown  in  Figure  57.  The  temperatures  appear  very  similar  in  that  there  is  only  a 
narrow  region  that  is  subjected  to  a  higher  temperature.  The  Szmerekovsky  model  tends 
to  more  uniformly  heat  the  material  and  also  to  have  a  higher  temperature.  This  is  due  to 
the  smaller  mesh  utilized  in  the  Szmerekovsky  model,  and  the  higher  energy  present  in 
the  Szmerekovsky  model,  that  must  be  turned  into  heat. 
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Figure  54  -  Velocity  Profile  -  Szmerekovsky 
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Figure  55  -  Velocity  Profile  -  Rod  Impact  Model 
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Figure  56  -  Temperature  Profile  -  Szmerekovsky 
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Figure  57  -  Temperature  Profile  -  Rod  Impact  Model 


97 


Additionally,  tracer  points  were  used  to  investigate  the  exact  response  of  specific 
points  of  material  throughout  the  deformation  event.  Points  were  chosen  at  a  distance  of 
0.01  cm  above  the  material  interface,  so  that  a  good  representation  of  the  boundary  layer 
in  the  rod  could  be  seen.  Since  gouges  were  seen  to  occur  at  x-positions  between  4.4  cm 
and  4.8  cm,  fifty  tracer  points  were  spaced  between  4.5  cm  and  5.0  cm.  A  diagram  of  the 
tracer  points  is  provided  in  Figure  58. 


x=4.5  x=5-° 

Target 

Figure  58  -  Tracer  Point  Locations  ( Rod  Impact  Model) 

As  deformation  occurred,  temperature  and  pressure  were  recorded,  to  investigate 

how  these  two  variables  vary  as  gouging  occurs.  Three  tracer  points  were  chosen  to 

study,  based  upon  their  movement  during  the  gouging  event.  The  tracer  points  chosen 

were  point  number  one,  which  is  located  at  x=4.5  cm  and  y=0.01  cm,  point  number  10, 

which  is  located  at  x=4.59  cm  and  y=0.01  cm,  and  point  number  20,  which  is  located  at 

x=4.69  cm  and  y=0.01  cm.  The  reason  these  points  were  chosen  will  be  seen  shortly. 

After  5.0  ps,  the  rod  impact  model  has  deformed  in  a  manner  more  reminiscent 
of  wear  than  gouging,  as  can  be  seen  in  Figure  59.  However,  to  the  far  right,  the  hump 
described  earlier  is  beginning  to  form  and  grow.  At  this  point  in  the  impact  event,  the 
tracer  points  have  not  been  moved  very  much  from  their  initial  positions,  with  the 
exception  of  the  tracer  points  around  tracer  point  number  27. 
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Gouging  has  fully  started  by  5.25  jlxs,  as  seen  in  Figure  60.  The  hump  has 
increased  in  size,  and  the  rod  material,  which  is  moving  to  the  right,  has  begun  to  overrun 
the  hump,  fonning  the  classic  jet  of  material.  The  hump  has  increased  in  size  to  0.01  cm 
in  height,  and  has  started  to  displace  tracer  point  20  upwards,  as  the  hump  continues  to 
move  to  the  left.  Points  1  and  10  have  not  been  displaced  visibly  at  this  time. 

At  5.25  ps,  the  hump  of  material  caused  by  the  impact  has  begun  to  form  the 
classic  material  jet  seen  in  previous  work.  By  5.5  ps,  jetting  is  clearly  the  dominant 
feature  of  deformation,  see  Figure  61.  At  this  time,  material  from  the  rod  is  clearly  being 
pulled  down  into  the  target  material,  as  shown  by  the  dip  in  the  tracer  points.  It  appears 
that  point  10  is  being  dragged  downward  towards  the  lower  jet  of  material.  Point  20 
however  has  not  been  involved  in  the  lower  jet,  and  appears  to  be  passing  above  the 
upper  jet  of  material. 

The  features  which  were  beginning  to  be  seen  at  5.5  ps,  continue  to  develop  at 
5.75  ps,  as  shown  in  Figure  62.  Tracer  point  10  is  clearly  being  pulled  down  into  the  jet 
of  VascoMax  300  material,  which  is  impinging  into  the  1080  steel  target.  Additionally,  it 
appears  that  point  1  may  also  be  dragged  downwards  into  the  lower  jet.  Point  20  has 
continued  to  pass  above  the  upper  jet  and  is  nearly  out  of  the  deformation  area. 

By  6.05  ps,  the  gouging  event  has  nearly  dragged  tracer  point  1  down  into  the 
lower  material  jet.  Point  10  appears  to  be  at  about  0.01  cm  above  the  original  interface, 
as  compared  to  0.003  cm  as  in  Figure  62.  This  shows  that  the  event  is  drawing  to  a  close, 
because  the  rod  material  is  beginning  to  ride  up  out  of  the  target  material,  which 
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eventually  would  lead  to  the  materials  separating  and  the  gouging  event  ending.  Tracer 
point  20  has  passed  beyond  the  gouge. 
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Figure  59  -  Tracer  Points  at  5.0  /us 
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Figure  60  -  Tracer  Points  at  5.25  /us 
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Figure  61  -  Tracer  Points  at  5.5  /us 
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Figure  62  -  Tracer  Points  at  5. 75  /us 
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Figure  63  -  Tracer  Points  at  6.05  /us 


The  previous  five  figures  present  a  physical  basis  for  understanding  what  is 
happening  at  each  tracer  point  at  a  given  time  step  during  deformation.  From  these 
figures,  temperature  and  pressure  will  be  looked  at  and  discussed. 

Tracer  point  1  was  the  tracer  point  farthest  from  the  initial  point  of  impact,  which 
was  at  4.8  cm.  Since  this  point  was  farther  from  the  impact  region  than  the  others,  it 
would  be  assumed  that  it  would  take  longer  for  the  stress  wave  to  reach  this  point  and 
produce  any  effects.  This  is  shown  to  be  accurate  below  about  5.7  ps,  since  the  values  of 
pressure  and  temperature  seen  in  Figure  64  are  nearly  constant.  Referring  to  Figure  62,  it 
is  seen  that  around  5.7  ps,  tracer  point  1  begins  to  descend  into  the  lower  material  jet.  At 
about  6.0  ps,  the  pressure  at  point  1  spikes  to  nearly  25.0  GPa,  at  this  point,  the  stress 
wave  affecting  point  1  is  reflecting  around  through  the  material  jet,  which  explains  the 
somewhat  erratic  nature  of  the  pressure  at  this  point.  At  a  peak  pressure  of  25.0  GPa,  the 
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VascoMax  300  steel  is  being  plastically  deformed,  since  the  yield  strength  of  the  steel  is 
only  14.47  GPa.  The  temperature  spikes  at  the  same  point  in  time  as  the  pressure,  and 
reaches  a  maximum  temperature  of  about  1300  K.  At  these  temperatures,  Szmerekovsky 
predicted  that  the  steel  would  form  austenite,  which  was  found  to  occur  in  recovered 
gouges  from  the  HHSTT.[2]  Temperatures  due  to  the  impact  were  not  high  enough  to 
melt  the  VascoMax  300  steel,  which  has  a  melting  point  of  23 10  K.[l] 
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Figure  64  -  Tracer  Point  1  (a)  Pressure  (b)  Temperature 
Tracer  point  10  reacted  to  the  impact  simulation  in  much  the  same  way  that  tracer 

point  1  did.  Tracer  point  10  showed  an  earlier  reaction  to  the  impact,  as  is  seen  in  Figure 

65.  At  about  5.5  ps,  the  pressure  and  temperature  both  begin  a  rapid  rise  to  their  peak 

values.  Referring  back  to  Figure  61,  it  is  seen  that  about  this  time  is  when  tracer  point  10 

begins  to  move  down  into  the  lower  material  jet.  At  point  10,  the  pressure  approaches 

33.0  GPa,  while  the  temperature  approaches  1300  K.  This  temperature  also  means  that 

tracer  point  10  is  austenite  steel.  The  peak  pressure  and  temperature  occurs  at  6.0  ps, 

which  corresponds  in  Figure  63  to  the  time  when  tracer  point  10  is  moving  up  the  upper 

material  jet.  At  this  point  in  the  gouging  event,  the  rod  is  still  moving  to  the  right  at  4.8 
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km/s.  With  tracer  point  10  in  the  position  that  it  is  in,  it  will  be  under  a  tremendous 
nonnal  load  due  to  the  velocity  of  the  rod,  thus  causing  the  very  high  pressure. 
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Figure  65  -  Tracer  Point  10  (a)  Pressure  (b)  Temperature 
The  last  tracer  point  investigated,  point  20,  was  seen  in  Figure  59  through  Figure 

63  to  pass  generally  above  the  gouge,  and  not  get  dragged  down  into  the  lower  material 

jet.  As  can  be  seen  in  Figure  66,  this  fact  leads  to  some  significantly  different  pressure 

and  temperature  histories.  The  pressure  and  temperature  both  begin  their  rapid  rise  at 

around  5.1  ps,  which  corresponds  roughly  to  the  point  at  which  the  hump  begins  to 

deflect  tracer  point  20  up.  The  peak  pressure  of  approximately  17  GPa  results  in  much 

less  plastic  deformation  in  this  area,  which  is  why  gouging  did  not  begin  earlier  in  the 

event  than  it  did.  The  rapid  drop  off  in  pressure  also  shows  that  all  of  the  plastic 

deformation  at  point  20  occurs  over  a  very  short  time  period.  The  temperature  graph 

leads  to  much  the  same  conclusions,  due  to  the  plateau  starting  around  5.5  ps.  From 

Figure  61,  it  was  seen  that  by  5.5  ps,  tracer  point  20  was  beyond  the  crest  of  the  material 

jet,  and  therefore  nearly  done  with  the  gouging  event.  Once  it  was  beyond  the  crest  of 
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the  jet  it  should  stay  a  fairly  constant  temperature,  because  there  is  no  more  directly 
applied  stress,  only  stress  due  to  reflections. 

L  POINT  20 


0  12  3  4  5  6  7 

TIME  Us) 

4.8:<.6em  circular  ftp  —  4.8km/s  by  3.2m/s 

(b) 

Figure  66  -  Tracer  Point  20  (a)  Pressure  (b)  Temperature 

In  addition  to  the  simulation  of  the  1500  m/s  model  developed  by  Szmerekovsky, 
a  3000  m/s  model  was  also  simulated.  In  this  case,  the  simulated  velocity  of  the  rod 
impact  specimen  was  9617  m/s  horizontal  by  6.4  m/s  vertical.  This  model  also  showed 
the  appearance  of  wear,  however  gouging  was  never  seen  to  occur.  The  projectile 
deformed  as  seen  in  Figure  67  where  the  wear  can  be  seen  all  along  the  rod  material  on 
top  and  the  rail  material  on  the  bottom.  The  mesh  utilized  in  this  simulation  was  0.0050 
cm  on  a  side. 

The  rod  impact  model  fired  at  9.6  km/s  did  not  develop  gouging.  This  model  did 
develop  some  of  the  characteristics  that  were  seen  in  the  slower  impact  model,  but 
gouging  never  did  occur  as  in  the  4.8  km/s  model.  The  deviatoric  stress  developed  a  high 
stress  region,  in  much  the  same  way  as  in  the  slower  model,  Figure  51,  however  gouging 
was  not  initiated  in  the  9.6  km/s  model,  as  seen  in  Figure  68.  The  materials  were  never 
able  to  gouge,  because  the  velocity  of  impact  was  much  greater  than  the  material  sound 
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speed,  and  the  stress  wave  was  never  able  to  propagate  ahead  of  the  material  which  is 
moving  at  9.6  km/s. 
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Figure  67  -  Rod  Impact  Model  at  9. 6  km/s 
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Figure  68  -  Deviatoric  Stress  -9.6  km/s  Model  -  No  Gouge 
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The  pressure  found  in  the  9.6  km/s  model  showed  a  similar  high-pressure  core  as 
seen  in  the  4.8  km/s  model.  The  gouge  in  the  4.8  km/s  model  was  initiated  when  a  high- 
pressure  core  was  formed,  and  rapid  plastic  deformation  began  to  occur,  as  seen  in  Figure 
69.  A  similar  high-pressure  core  was  also  seen  in  the  9.6  km/s  model,  with  pressures 
very  nearly  equal  to  those  seen  in  the  4.8  km/s  model.  This  high-pressure  core  for  the  9.6 
km/s  model  is  seen  in  Figure  70.  The  only  difference  is  that  in  the  slower  model,  the 


high-pressure  core  causes  the  gouging  to  occur,  while  the  faster  model  it  does  not. 
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Figure  69  -High  Pressure  Core  at  Gouge  Initiation  (4.8  km/s  Model) 
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Figure  70  -  High  Pressure  Core  (9.6  km/s  Model)  -  No  Gouge 


4.3  Experimentally  Verifiable  Rod  Impact  Model 

The  problem  with  the  two  models  discussed  earlier,  the  1500  m/s  and  3000  m/s 
models,  is  that  to  match  Szmerekovsky’s  models  of  the  impact,  it  is  necessary  to  launch 
the  projectiles  at  4.8  km/s  and  9.6  km/s,  respectively.  Facilities  available  on  Wright 
Patterson  AFB,  capable  of  shooting  impact  specimens,  currently  provide  the  ability  to 
reach  approximately  2133  m/s,  which  is  much  lower  than  is  needed  to  test  the  models 
developed  here.  This  led  to  the  design  of  another  model,  based  upon  the  experimental 
capabilities  available. 
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Simply  performing  a  simulation  using  the  velocity  of  the  gun  available  does  not 
provide  much  insight  into  how  the  model  relates  to  the  HHSTT.  Therefore,  it  was 
decided  to  use  the  Buckingham  Pi  process  in  the  reverse  manner  to  determine  what  the 
rod  impact  model  at  2133  m/s  represents  in  terms  of  the  HHSTT.  To  begin,  the  size  of 
the  projectile  was  kept  the  same  as  for  both  of  the  previous  models.  By  keeping  the 
diameter  equal  to  0.6  cm  and  the  length  equal  to  4.8  cm,  scaling  was  maintained  with 
both  the  Szmerekovsky  model  and  the  original  HHSTT  model.  To  maintain  the  velocity 
scaling,  the  original  n3  was  used,  which  was  0. 125.  The  horizontal  velocity  was  then 
assumed  to  be  2133  m/s,  and  from  this,  the  Pi  parameters  used  in  scaling  down  the 
Szmerekovsky  model  were  calculated  for  this  new  model.  These  values  are  seen  in  Table 
7. 


Table  7  -  Experimentally  Verifiable  Rod  Impact  Model  Pi  Parameters 


Experimentally 
Reproducible  Rod 
Impact  Model 

Ux 

2133  m/s 

Uy 

-2  m/s 

n. 

0.125 

n3 

-6.667x1  O'4 

n4 

2.221 

n5 

0.500 

n6 

0 

n7 

1.432 

n8 

180.5 

n9 

71.066 

The  goal  was  to  convert  these  values  back  into  an  equivalent  experiment  at  the 
HHSTT.  This  was  accomplished  by  going  back  to  the  definition  of  the  Pi  parameters 
from  Section  3.3.  From  the  original  development,  EE  was  shown  to  be, 
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Similar  equations  exist  for  the  elastic  modulus  (n8)  and  the  shear  modulus  (rig),  however 
they  do  not  add  anything  to  the  solution,  because  the  only  differences  between  the  three  is 
the  first  term  on  the  right,  which  we  assume  are  constants  with  the  material.  The  second 
tenn  on  the  right  of  the  three  definitions  is  exactly  the  same;  therefore  it  is  pointless  to 
solve  the  second  two  equations.  The  mass,  m,  in  the  above  equation  is  the  combined 
mass  of  the  equivalent  slipper  and  the  sled  mass.  From  this  equation,  a  few  assumptions 
are  required  as  to  the  conditions  at  the  HHSTT.  The  first  assumption  is  that  the  same 
slipper  is  in  use  for  the  equivalent  Holloman  model  as  is  used  for  the  original  HHSTT 
model  developed  by  Szmerekovsky.  This  means  that  the  yield  strength  and  mass  of  the 
equivalent  slipper  are  known,  as  well  as  the  length  in  the  numerator.  This  leaves  the 
mass  of  the  equivalent  sled  and  the  equivalent  velocity  as  the  only  two  unknowns.  For 
this  research,  the  worst  case  of  a  mass  of  800  kg  was  assumed,  and  therefore  the 
equivalent  velocity  was  found  to  be, 
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From  this  last  equation,  and  with  the  assumptions  presented  in  the  previous  discussion, 
the  equivalent  HHSTT  experiment  can  be  shown  to  have  the  parameters  seen  in  Table  8, 
and  a  comparison  between  the  rod  impact  model  and  the  equivalent  HHSTT  experiment 
are  shown  pictorially  in  Figure  7 1 . 
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Table  8  -  Equivalent  HHSTT Experiment 


Rod 

Impact 

Model 

Equivalent 

HHSTT 

Experiment 

ux  (m/s) 

2133 

669.0 

uv(m/s) 

-1.4 

-0.446 

Length  (cm) 

4.8 

20.32 

Diameter  (cm) 

0.6 

2.54 

Sled  Mass  (kg) 

Not 

Applicable 

800 

Experimentally 

Verifiable  Rod  Impact  Equivalent  HHSTT 

Model 


Figure  71  -  Equivalent  HHSTT  Experimental  Sled 


Simulation  using  these  parameters  predicted  deformation  as  seen  in  Figure  72. 
Deformation  seen  is  more  reminiscent  of  wear  as  described  earlier. 

Simulation  of  an  800  kg  sled  at  669  m/s  did  not  show  the  characteristics  of 
gouging.  It  was  decided  to  alter  the  rod  impact  model  parameters  to  increase  the 
equivalent  velocity  of  impact  that  the  HHSTT  would  see.  The  two  parameters  that  can  be 
changed  are  the  diameter  of  the  impact  rod,  and  the  mass  assumed  for  the  sled.  A  variety 
of  combinations  were  investigated  by  changing  the  rod  diameter  and  the  sled  mass. 

Table  9  shows  a  few  of  the  options  that  were  developed. 
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Figure  72  -  Rod  Impact  Model  Simulation  at  2.13  km/s  -  No  Gouge 


Table  9  -  Equivalent  HHSTT  Velocity  by  Varying  Input  Parameters 


Option 

Rod  Diameter 
(cm) 

Sled  Mass 
(kg) 

Equivalent  HHSTT 
Velocity  (m/s) 

1 

.6 

600 

769.7 

2 

.6 

300 

1073.0 

3 

.5 

600 

843.1 

4 

.5 

300 

1175.0 

5 

.4 

600 

942.7 

6 

.4 

300 

1314.0 

7 

.3 

600 

1089.0 

8 

.3 

300 

1517.0 

It  was  decided  to  run  only  the  last  option  which  has  the  highest  equivalent 
velocity  at  the  HHSTT,  because  it  has  been  witnessed  in  the  field  that  generally  gouging 
does  not  occur  below  1500  m/s.  The  equivalent  HHSTT  test  sled  is  shown  in  Figure  73. 
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The  results  for  this  simulation  are  shown  in  Figure  74.  Again,  at  this  velocity,  it  appears 
that  only  wear  occurs,  and  a  gouge  does  not  form. 
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Figure  73  -  Equivalent  HHSTT Experimental  Sled  (300  kg) 
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Figure  74  -  Rod  Impact  Model  at  2133  m/s  -  0.3  cm  diameter  -  Simulates  300  kg  Sled 
Simulations  involving  various  test  sled  masses  and  velocities  have  shown  that 

there  is  no  single  value  of  sled  mass  or  velocity  that  can  guarantee  gouging  would  occur. 

It  was  shown  that  a  sled  at  1500  m/s  with  a  mass  of  800  kg  gouges,  while  a  sled  at  the 

same  velocity,  but  a  mass  of  300  kg  will  not  gouge.  Similarly,  a  sled  with  a  mass  of  800 
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kg  will  gouge  at  1500  m/s,  but  will  not  gouge  at  669  m/s.  These  facts  lead  to  the 
assumption  that  there  is  a  relationship  between  velocity  and  mass  that  would  predict 
when  gouging  would  occur. 
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Chapter  5  -  Conclusions  and  Recommendations 


5.1  Conclusions 

The  goal  of  this  research  was  to  develop  a  ballistic  impact  model  that  would 
simulate  gouging  as  seen  at  the  Holloman  High  Speed  Test  Track.  The  model  developed 
must  be  experimentally  verifiable,  using  the  compressed  gas  guns  available  at  Wright 
Patterson  Air  Force  Base.  The  hydrocode  computer  program  CTH  was  used  to  simulate 
the  ballistic  impacts,  to  investigate  whether  the  gouging  phenomena  was  present.  To 
reach  the  main  goal  of  this  thesis,  a  series  of  different  models  were  used. 

The  first  model  utilized  was  a  simulation  of  a  Taylor  impact  test.  This  simulation 
was  performed  to  verity  the  parameters  used  in  future  CTH  simulations.  The  Taylor  test 
was  chosen  as  a  verification  tool,  because  of  its  simplicity,  and  the  level  of  understanding 
of  the  principles  behind  it.  Additionally,  the  data  and  pictures  obtained  from  the  article 
by  Jones  et  al.  [14]  were  perfectly  suited  to  verifying  CTH  simulations. 

The  baseline  CTH  model  was  based  upon  the  model  used  by  Szmerekovsky, 
making  use  of  the  options  that  he  used.  The  simulation  appeared  to  give  fairly  accurate 
deformations,  see  Figure  27(b).  The  options  that  were  altered  generally  gave  expected 
results.  It  was  imagined  that  changing  the  yield  strength  in  mixed  cells  to  zero  in  mixed 
cells  would  give  unrealistically  low  material  strength,  which  was  shown  to  be  the  case  by 
the  specimen  appearing  to  erode  away  at  the  interface  due  to  material  fracturing  from  the 
bulk  of  the  material  and  separating.  Changing  the  manner  in  which  materials  and 
pressures  in  mixed  cells  were  allocated  also  changed  the  solution,  except  in  this  case  it 
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was  for  the  better,  which  was  expected,  because  the  baseline  model  did  not  allow  for 
pressure  relaxation  within  a  cell,  which  is  needed  for  the  conservation  equations.  Giving 
voids  strength,  by  not  compressing  them  from  the  material  first,  created  a  situation  where 
fracture  occurred  unrealistically,  again  causing  erosion  that  was  not  seen  experimentally. 
The  last  option,  the  slide  line  approximated  the  impact  surface  as  being  frictionless,  and 
therefore  only  plasticity  played  a  part  at  the  interface.  This  plasticity  at  the  interface 
developed  deviatoric  stresses  farther  from  the  interface.  This  option  gave  good  results  as 
well,  and  matched  the  experimental  conditions  that  the  end  of  the  rod  and  the  target  were 
both  polished. 

After  verifying  that  the  correct  input  options  were  being  used,  the  rod  impact 
model  was  developed  and  sized  using  the  Buckingham  Pi  theorem.  The  rod  impact 
model  was  used  to  simulate  the  1500  m/s  HHSTT  experiment,  which  was  simulated 
realistically  by  Szmerekovsky.  This  model  was  shown  to  develop  both  wear  and  gouging 
of  the  materials.  The  wear  developed  prior  to  the  gouging,  which  is  accurate;  because 
gouging  is  actually  a  case  where  the  wear  builds  up  to  a  critical  point,  at  which  the 
materials  fonn  a  hump  and  jet  into  each  other. 

The  gouges  developed  using  the  rod  impact  model,  Figure  42,  and  the 
Szmerekovsky  model,  Figure  41,  do  not  appear  exactly  the  same,  due  in  large  part  to  the 
scale  of  each  simulation.  The  gouge  in  the  rod  impact  model  is  seen  to  be  approximately 
0.03  cm  tall,  while  that  in  the  Szmerekovsky  model  is  0. 1  cm  tall.  Within  the  interaction 
zone,  the  mesh  of  the  rod  impact  model  is  0.0050  cm,  while  the  Szmerekovsky  model 
uses  a  mesh  of  0.0025  cm.  In  effect,  there  are  only  six  cells  vertically  in  the  gouged  zone 
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of  the  rod  impact  model,  while  the  Szmerekovsky  model  has  forty  cells  in  the  gouged 
area.  This  explains  the  lack  of  smoothness  in  the  rod  impact  model,  which  is  seen  in  the 
Szmerekovsky  model. 

The  rod  impact  model  did  a  very  good  job  of  modeling  the  gross  defonnation  seen 
in  the  formation  of  a  gouge,  Figure  42.  Exact  dimensionality  was  not  achieved  with  the 
rod  impact  model,  due  to  the  differences  in  scale  between  it  and  the  Szmerekovsky 
model.  The  rod  impact  model  also  represented  the  high-pressure  core  seen  by 
Szmerekovsky  very  well.  The  very  high  shear  stress  gradient  on  the  crest  of  the  jet  seen 
in  Szmerekovsky ’s  work,  was  simulated  quite  well  by  the  rod  impact  model. 
Szmerekovsky  showed  that  the  material  within  and  around  a  gouge  tends  to  flow  up  and 
along  the  jetting  material.  This  feature  of  material  flow  was  also  seen  using  the  rod 
impact  model.  Lastly,  the  rod  impact  model  matched  the  plastic  defonnation  seen  in  the 
Szmerekovsky  model  quite  well. 

The  time  scale  developed  through  the  Buckingham  Pi  relationship  did  not  match 
very  well  between  the  two  models.  It  appeared  that  gouging  occurred  at  the  same  time 
for  both  the  Szmerekovsky  model  and  the  rod  impact  model,  whereas  it  should  have 
occurred  at  0.44  ps  in  the  rod  impact  model,  as  opposed  to  the  6.0  ps  that  it  actually 
occurred.  There  are  many  factors  that  could  cause  this  parameter  to  be  off.  First,  the 
time  of  interest  may  not  actually  be  taken  from  the  beginning  of  the  simulation,  it  may  be 
the  taken  from  the  initiation  of  gouging.  The  mass  of  material  directly  above  the  area 
that  gouges  may  also  play  a  role,  in  adding  inertia  directly  into  the  gouge.  The 
Szmerekovsky  model  is  much  thicker  than  the  rod  impact  model,  43  cm  for 
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Szmerekovsky  and  0.6  cm  for  the  rod  impact  model.  All  of  the  possibilities  go  back  to 
the  fact  that  only  ten  invariants  were  found  in  the  Buckingham  Pi  analysis  of  the  two 
models.  If  one  were  able  to  put  all  possible  physical  variables  between  the  two  models 
into  a  single  Buckingham  Pi  analysis,  then  the  rod  impact  model  developed  would 
exactly  match  the  Szmerekovsky  model.  Overall  though,  the  rod  impact  model  develops 
an  excellent  representation  of  gouging,  even  if  time  scales  are  not  exact. 

Gouging  never  occurred  in  the  simulation  of  the  3000  m/s  HHSTT  sled  run.  This 
is  due  to  the  fact  that  as  the  materials  interact,  the  relative  velocity  between  them  is  9600 
m/s,  while  the  velocity  at  which  energy  can  be  dissipated  (the  material  sound  speed),  is 
only  4930  m/s.  This  means  that  the  materials  are  interacting  at  almost  twice  the  speed  at 
which  any  energy  can  be  transferred.  Because  of  this,  the  materials  will  instantaneously 
fracture  due  to  the  increase  of  energy  in  the  cells  at  the  interface.  This  disintegration  of 
material  has  been  called  erosion.  Due  to  this  fact,  the  rod  impact  model  developed  here  is 
incapable  of  modeling  what  occurs  at  the  HHSTT  at  a  velocity  of  3000  m/s.  It  may  be 
possible  to  artificially  develop  gouging  at  this  velocity  if  yield  strength  is  adjusted  higher. 

The  last  two  models  were  designed  to  simulate  conditions  that  are  actually 
possible  to  duplicate  in  laboratory  facilities  at  Wright  Patterson  AFB.  The  first  model 
simulated  an  actual  sled  arrangement  utilized  at  the  HHSTT,  but  the  equivalent  HHSTT 
velocity  was  found  based  upon  the  velocity  possible  at  Wright  Patterson  AFB.  Using  the 
Buckingham  Pi  theorem,  along  with  the  invariants  developed  in,  the  equivalent  HHSTT 
velocity  was  found  to  be  669.0  m/s.  This  velocity  is  much  lower  than  the  velocity  where 
gouging  is  generally  found,  so  it  was  expected  that  a  gouge  would  not  form.  Simulations 
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proved  that  at  this  velocity  and  model  size,  a  gouge  did  not  form,  however  the 
characteristics  associated  with  wear  were  present. 

At  velocities  lower  than  approximately  1500  m/s,  the  HHSTT  does  not  generally 
see  gouging.  Deformations  that  are  seen  in  this  velocity  range  tend  to  be  described  as 
wear,  which  was  described  earlier  as  “progressive  damage  to  a  surface  caused  by  relative 
motion  with  respect  to  another  substance.’’[15]  Wear  was  witnessed  for  the  relative 
velocity  of  669.0  m/s,  as  modeled  using  a  rod  impact  velocity  of  2 133  m/s.  The  likely 
reason  that  gouging  did  not  occur  at  this  velocity  is  that  the  kinetic  energy  of  a  sled 
traveling  at  1500  m/s  is  five  times  greater  than  that  of  a  sled  traveling  669.0  m/s.  It  was 
clear  that  to  initiate  gouging,  a  parameter  or  parameters  would  have  to  be  changed,  to 
increase  the  relative  velocity  into  the  range  where  gouging  occurs. 

A  combination  of  decreasing  the  assumed  weight  of  the  test  sled,  and  also 
decreasing  the  diameter  of  the  impacting  rod  increased  the  relative  velocity  of  the 
simulated  HHSTT  sled  up  to  1517  m/s.  Simulation  using  these  parameters  did  not 
develop  gouging  either.  This  leads  to  the  conclusion  that  velocity  is  not  the  only  factor 
important  to  the  gouging  phenomena.  The  other  factor  that  is  likely  to  cause  gouging  is 
the  mass  of  the  sled.  In  order  to  reach  a  relative  velocity  of  1500  m/s,  it  was  necessary  to 
assume  that  the  test  sled  had  a  mass  of  only  300  kg.  A  sled  of  mass  800  kg  would  have 
2.67  times  the  kinetic  energy  of  a  sled  with  a  mass  of  300  kg.  The  higher  kinetic  energy 
of  the  larger  sled  would  deposit  more  energy  into  the  rail  than  would  the  smaller  sled, 
causing  higher  temperatures  and  increased  plasticity. 
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The  rod  impact  model  was  shown  to  simulate  gouging  exceptionally  well  in  a 
global  sense.  Overall,  the  defonnation  seen  with  the  rod  impact  model  shows  that  a 
gouge  is  fully  developed  using  the  simplified  model.  Additionally,  some  of  the  features 
seen  in  the  pressure  contours  and  shear  stress  contours  showed  excellent  correlation. 
Differences  between  the  Szmerekovsky  model  and  the  rod  impact  model  become 
apparent  when  values  of  pressure,  stress,  temperature,  etc.  are  investigated.  Numerical 
differences  between  the  models  suggest  that  perhaps  in  the  Buckingham  Pi  process, 
important  dimensioned  quantities  were  left  out.  In  scaling  the  rod  impact  model,  13 
dimensioned  quantities  were  used.  By  using  more  dimensioned  quantities,  perhaps  better 
numerical  agreement  could  be  achieved. 

In  conclusion,  the  model  developed  in  this  research  does  simulate  the  formation 
of  gouging  in  a  high  velocity  oblique  impact  experiment.  The  gouge  developed  is  very 
good  in  a  global  sense,  meaning  that  gross  deformation  is  modeled  well,  however  local 
values  of  state  variables  are  not  exact.  This  leads  to  the  conclusion  that  the  model 
developed  here  is  useful  for  investigating  gouging,  but  at  the  moment  is  not  useful  in 
design  of  slipper/rail  arrangements.  The  equivalent  HHSTT  models  that  can  be 
investigated  using  facilities  available  at  Wright  Patterson  AFB,  were  shown  to  develop 
characteristics  of  wear,  which  means  that  the  equivalent  model  could  be  useful  in  better 
understanding  how  wear  affects  the  HHSTT. 
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5.2  Recommendations  for  Future  Research 


It  was  assumed  by  Szmerekovsky  that  the  mass  of  the  sled  was  evenly  distributed 
between  the  four  slippers.  This  may  be  fairly  precise  when  the  test  sled  is  at  rest,  but  it 
may  not  hold  any  validity  during  an  actual  experiment.  A  combination  of  aerodynamic 
forces,  rail  characteristics,  and  sled  dynamics  can  lead  to  the  sled  moving  with  a  yaw,  a 
pitch,  or  a  roll.  Any  of  these  movements  could  alter  the  weight  distribution  between  the 
four  slippers,  possibly  causing  an  impact  with  the  rail  where  only  one  slipper  was  in 
contact,  thereby  applying  the  entire  load  to  a  single  slipper,  increasing  four  fold  the 
kinetic  energy  which  was  used  in  both  the  Szmerekovsky  model  and  the  rod  impact 
model.  Therefore,  the  validity  of  the  manner  in  which  the  sled  mass  is  divided  should  be 
investigated. 

The  ultimate  goal  of  this  research  was  to  develop  a  simplified  model  of  gouging, 
which  could  be  experimentally  verified  to  be  correct.  This  goal  was  accomplished,  in 
that  a  model  was  developed  that  showed  the  development  of  a  gouge  similar  to  the  gouge 
seen  in  research  by  Szmerekovsky  and  others.  Experimental  verification  of  the  rod 
impact  model  with  a  diameter  of  0.6  cm  and  a  length  of  4.8  cm,  fired  at  a  velocity  of  4.8 
km/s  would  provide  verification  that  the  model  developed  herein  is  accurate.  Achieving 
this  velocity  would  require  use  of  a  different  compressed  air  gun  than  what  is  available 
for  use  at  Wright  Patterson  AFB.  To  further  verify  the  rod  impact  model,  it  would  be 
necessary  to  develop  a  model  that  matches  time  scales  with  the  Szmerekovsky  model. 
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Additional  simulations  using  a  variety  of  velocities  and  masses  would  be  very 
valuable  in  further  verifying  what  combinations  of  mass  and  velocity  develop  gouging. 
Further  simulations  would  be  able  to  better  clarify  where  the  lines  between  no  gouging 
and  gouging  are  located.  This  graph  could  then  be  compared  with  observations  made  by 
the  HHSTT  engineers  of  when  gouging  occurs  in  physical  experiments.  If  possible, 
gouges  seen  in  simulation  should  be  experimentally  verified. 

In  general,  the  test  track  at  Holloman  AFB  has  an  epoxy  coating  on  the  rails, 
which  serves  as  a  type  of  lubricant  to  prevent  gouging.  Addition  of  a  coating  to  the  target 
material  of  the  models  developed  here  would  allow  characterization  of  the  manner  in 
which  coatings  would  affect  the  rod  impact  model.  These  simulations  could  be  compared 
to  the  work  by  Szmerekovsky  that  studied  the  affects  of  a  coating  on  the  gouging 
phenomena. 
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Appendix  1 


CTH  Input  Deck  -  Taylor  Impact  Model 


*eor*  genin 

Taylor  Test:  Copper  on  4340  Steel,  V=176  m/s 
control 
inmp 
ep 

vpsave 

endcontrol 

mesh  *  Define  the  mesh  using  this  area 

block  1  geom=2dc  type=e 
x0=0.0 

xl  n=80  w=2  dxf=0.025 
x2  n=40  w=4dxf=0.1 
endx 


y0=-20.0 

yl  n=150  w=15  dyf=0.1 
y2  n=480  w=12  dyf=0.025 
endy 
endb 
endmesh 

insertion  of  material 
block  1 

package  topblock 
material  1 
numsub  50 
yvel  -176e2 
insert  box 
pi  0.0  0.01 
p2  0.381  5.725 
endinsert 
endpackage 
package  bottomblock 
material  2 
numsub  50 
insert  box 


*  Define  the  areas  within  mesh  that 

*  contain  material,  apply  velocities 


pi  0.0,  0.0 
p2  6, -17.5 
endinsert 
endpackage 
endblock 
endinsertion 
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edit 
block  1 
expanded 
endblock 
endedit 
tracer 

add  0.381,  0.01  to  0.381,  5.725  number=60 
endt 
eos 

MAT1  SES  COPPER 
MAT2  SES  STEEL  4340 
endeos 
epdata 
mix  3 

matep  1  jo=l  *  copper 
matep  2  jo=4340_TEMP_MART 
slide  1  2 
vpsave 
lstrain 
endep 


*  This  defines  tracer  points 


*  Set  equations  of  state  for 

*  materials  inserted  earlier 


*  Define  elastic-plastic 

*  material  properties 


*eor*  cthin  *  This  section  describes  how 

Taylor  Test:  Copper  on  4340  Steel,  V=176  m/s  *  CTH  will  run  problem 
control 
mmp 

tstop  =  100.0e-6  *  Defines  time  for  sim  to  stop 

endc 
Convct 
convection=l 
interface=high_resolution 
endc 
fracts 

pfracl  -10.0e9 
pfimix  -12.0e9 
pfvoid  -12.0e9 
endf 
edit 
shortt 

tim  0.0,  dt  =  1.0 
ends 
longt 

tim  0.0,  dt  =  1.0 
endl 
plott 

tim  0.0  dt  =  1.0e-6  *  Defines  when  to  record 


124 


*  data 


endp 

histt 

tim  0.0,  dt  =  1.0e-6 
htracer  all 
endh 
ende 

boundary 
bhy 
bl  1 

bxb  =  0  ,  bxt  =  1 
byb  =  1  ,  byt  =  1 
endb 
endh 
endb 

cellthermo 

mmp3 

ntbad  100000 
ende 


*  Sets  boundary  conditions 
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Appendix  2 


CTH  Input  Deck  -  Oblique  Rod  Impact  Model 

*eor*  genin 

V300/1080  4.8x.6cm  circular  tip  -  ~4.8km/s  by  3.2m/s 
* 

* 

* 

control 

mmp 

ep 

vpsave 

endcontrol 

* 

* 

* 

mesh 

block  1  geom=2dr  type=e  *  define  mesh 

x0=0.0 

xl  n=40  w=4.000  dxf=0.1000 
x2  n=60  w=0.600  dxf=0.0100 
x3  n=160  w=0.800  dxf=0.0050 
x4  n=60  w=0.600  dxf=0.0100 
x5  n=20  w=2.000  dxf=0.1000 
endx 

* 

y0=-4.0 

yl  n=20  w=2.000  dyf=0.1000 
y2  n=140  w=1.400  dyf=0.0100 
y3n=240w=  1.200  dyf=0.0050 
y4  n=40  w=0.400  dyf=0.1000 
endy 
endblock 

endmesh 

* 

* 

* 

insertion  of  material 

block  1  *  insert  cylindrical  rod  and  rotate  the  rod 

package  rod 
material  1 
numsub  100 
xvel  4808. 5e2 
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yvel  -3.21e2 
insert  box 
pi  0.0  0.0 
p2  4.8  0.6 
ppoint  4.8  0.3 
angle  -0.0382 
endinsert 
delete  circle 
center  4.8  0.3 
radius  0.3 
enddelete 
endpackage 

package  tip  *  insert  circular  tip  of  rod 

material  1 
numsub  100 
xvel  4808. 5e2 
yvel  -3.21e2 
insert  circle 
center  4.8  0.3 
radius  0.3 
endinsert 
endpackage 

package  target  *  insert  “rail”  material 

material  2 
numsub  100 
insert  box 
pi  0.0 -4.0 

p2  8.0  0.0 

endinsert 

endpackage 

endblock 

endinsertion 

* 

* 

* 

edit 

blockl 

expanded 

endblock 

endedit 

* 

* 

* 
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*  define  material  EOS 


eos 

MAT1  SES  STEELV300 
MAT2  SES  IRON 

endeos 
* 

* 

* 

epdata  *  define  constitutive  equations  of  materials 

mix  3 
matep  1 
st=19 

yield=14.47e9 
poisson=0.27 
matep  2 

johnson-cook=IRON 
yield=7.0e9 
poisson=0.28 
slide  1  2 
vpsave 
lstrain 

endep 

* 

* 

* 

* 

*eor*  cthin 
* 

V300/1080  4.8x.6cm  circular  tip  -  ~4.8km/s  by  3.2m/s 
* 

* 

* 

control 

mmp 

tstop=10.0e-6  *  define  stop  time  of  sim 

endcontrol 

* 

* 

* 

convct 

convection=l 

interface=high_resolution 

endconvct 

* 

* 

* 
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fracts 

pfracl=-10.0e9 

pfrac2=-2.5e9 

pfmix=-12.0e9 

pfvoid=-12.0e9 

endfracts 

* 

* 

* 

edit 

shortt 

tim  0.0,  dt=1.0 
ends 
longt 

tim  0.0,  dt=1.0 
endl 
plott 

tim  0.0,  dt=0.050e-6 
endp 

ende 

* 

* 

* 

boundary 
bhy 
bl  1 

bxb=l  ,  bxt=l 
byb=l  ,  bxt=l 
endb 
endh 

endb 

* 

* 

* 

vadd 

block=l 

tadd=0.0 

xvel=-4808.5e2  *  apply  a  velocity  to  keep  gouge  in  mesh 

endvadd 

* 

* 

* 

cellthermo 

mmp3 

ntbad  1000000 
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endc 

* 

* 

* 
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Pressure  History  Plots 


Appendix  3 
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All  points  beyond  point  40  are  constant  pressure  because  they  are  not  in  the  rod. 
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Temperature  History  Plots 


Appendix  4 
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All  points  beyond  point  40  are  constant  temperature  because  they  are  not  in  the  rod. 
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